#Loading required libraries

library(phyloseq)
# Define packages
## CRAN repositories
cran_packages <- c("bookdown", "knitr", "tidyverse", "plyr", "grid", "gridExtra", "kableExtra", "xtable", "ggpubr")
## Bioconductor repository
bioc_packages <- c("phyloseq", "dada2", "DECIPHER", "phangorn", "ggpubr", "BiocManager","DESeq2", "microbiome", "philr")
##  GitHub repository
git_source <- c("twbattaglia/btools", "gmteunisse/Fantaxtic", "MadsAlbertsen/ampvis2", "opisthokonta/tsnemicrobiota") # fuente/nombre del paquete
git_packages <- c("btools", "fantaxtic", "ampvis2", "tsnemicrobiota") # nombre del paquete

# Install CRAN packages
.inst <- cran_packages %in% installed.packages()
if(any(!.inst)) {
  install.packages(cran_packages[!.inst])
}
# Install Bioconductor packages
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
.inst <- bioc_packages %in% installed.packages()
if(any(!.inst)) {
  BiocManager::install(bioc_packages[!.inst])
}

## Bioconductor version 3.12 (BiocManager 1.30.16), R 4.0.3 (2020-10-10)

## Installing package(s) 'philr'

## also installing the dependency 'ggtree'

## Warning in .inet_warning(msg): installation of package 'ggtree' had non-zero
## exit status

## Warning in .inet_warning(msg): installation of package 'philr' had non-zero exit
## status

## Updating HTML index of packages in '.Library'

## Making 'packages.html' ... done
## Old packages: 'aplot', 'class', 'fansi', 'fftwtools', 'foreign', 'MASS',
##   'nlme', 'nnet', 'Rcpp', 'reticulate', 'spatial', 'tidytree'

# Install GitHub packages
.inst <- git_source %in% installed.packages()
if(any(!.inst)) {
  devtools::install_github(git_source[!.inst])
}

## Skipping install of 'btools' from a github remote, the SHA1 (fa21d4ca) has not changed since last install.
##   Use `force = TRUE` to force installation
## Skipping install of 'fantaxtic' from a github remote, the SHA1 (3342beee) has not changed since last install.
##   Use `force = TRUE` to force installation
## Skipping install of 'ampvis2' from a github remote, the SHA1 (f48b708f) has not changed since last install.
##   Use `force = TRUE` to force installation
## Skipping install of 'tsnemicrobiota' from a github remote, the SHA1 (cffbab72) has not changed since last install.
##   Use `force = TRUE` to force installation

# Loading required packages
sapply(c(cran_packages, bioc_packages, git_packages), require, character.only = TRUE)

## Loading required package: bookdown
## Loading required package: knitr
## Loading required package: tidyverse
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.6     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.4     ✓ stringr 1.4.0
## ✓ readr   2.1.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## Loading required package: plyr
## ------------------------------------------------------------------------------
## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)
## ------------------------------------------------------------------------------
## 
## Attaching package: 'plyr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## 
## The following object is masked from 'package:purrr':
## 
##     compact
## 
## Loading required package: grid
## Loading required package: gridExtra
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Loading required package: kableExtra
## 
## Attaching package: 'kableExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     group_rows
## 
## Loading required package: xtable
## Loading required package: ggpubr
## 
## Attaching package: 'ggpubr'
## 
## The following object is masked from 'package:plyr':
## 
##     mutate
## 
## Loading required package: dada2
## Loading required package: Rcpp
## Loading required package: DECIPHER
## Loading required package: Biostrings
## Loading required package: BiocGenerics
## Loading required package: parallel
## 
## Attaching package: 'BiocGenerics'
## 
## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB
## 
## The following object is masked from 'package:gridExtra':
## 
##     combine
## 
## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union
## 
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## 
## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, basename, cbind, colnames,
##     dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
##     grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
##     order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
##     rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
##     union, unique, unsplit, which.max, which.min
## 
## Loading required package: S4Vectors
## Loading required package: stats4
## 
## Attaching package: 'S4Vectors'
## 
## The following object is masked from 'package:plyr':
## 
##     rename
## 
## The following objects are masked from 'package:dplyr':
## 
##     first, rename
## 
## The following object is masked from 'package:tidyr':
## 
##     expand
## 
## The following object is masked from 'package:base':
## 
##     expand.grid
## 
## Loading required package: IRanges
## 
## Attaching package: 'IRanges'
## 
## The following object is masked from 'package:plyr':
## 
##     desc
## 
## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice
## 
## The following object is masked from 'package:purrr':
## 
##     reduce
## 
## The following object is masked from 'package:phyloseq':
## 
##     distance
## 
## Loading required package: XVector
## 
## Attaching package: 'XVector'
## 
## The following object is masked from 'package:plyr':
## 
##     compact
## 
## The following object is masked from 'package:purrr':
## 
##     compact
## 
## 
## Attaching package: 'Biostrings'
## 
## The following object is masked from 'package:base':
## 
##     strsplit
## 
## Loading required package: RSQLite
## Loading required package: phangorn

## Warning: package 'phangorn' was built under R version 4.0.5

## Loading required package: ape
## 
## Attaching package: 'ape'
## 
## The following object is masked from 'package:Biostrings':
## 
##     complement
## 
## The following object is masked from 'package:ggpubr':
## 
##     rotate
## 
## Loading required package: BiocManager
## Bioconductor version '3.12' is out-of-date; the current release version '3.14'
##   is available with R version '4.1'; see https://bioconductor.org/install
## Loading required package: DESeq2
## Loading required package: GenomicRanges
## Loading required package: GenomeInfoDb
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'matrixStats'
## 
## The following object is masked from 'package:plyr':
## 
##     count
## 
## The following object is masked from 'package:dplyr':
## 
##     count
## 
## 
## Attaching package: 'MatrixGenerics'
## 
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 
## Loading required package: Biobase
## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.
## 
## 
## Attaching package: 'Biobase'
## 
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## 
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## 
## The following object is masked from 'package:phyloseq':
## 
##     sampleNames
## 
## Loading required package: microbiome
## 
## microbiome R package (microbiome.github.com)
##     
## 
## 
##  Copyright (C) 2011-2020 Leo Lahti, 
##     Sudarshan Shetty et al. <microbiome.github.io>
## 
## 
## Attaching package: 'microbiome'
## 
## The following object is masked from 'package:SummarizedExperiment':
## 
##     coverage
## 
## The following object is masked from 'package:GenomicRanges':
## 
##     coverage
## 
## The following object is masked from 'package:phangorn':
## 
##     diversity
## 
## The following object is masked from 'package:Biostrings':
## 
##     coverage
## 
## The following objects are masked from 'package:IRanges':
## 
##     coverage, transform
## 
## The following object is masked from 'package:S4Vectors':
## 
##     transform
## 
## The following object is masked from 'package:ggplot2':
## 
##     alpha
## 
## The following object is masked from 'package:base':
## 
##     transform
## 
## Loading required package: philr

## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'philr'

## Loading required package: btools
## Loading required package: fantaxtic
## Loading required package: reshape2
## 
## Attaching package: 'reshape2'
## 
## The following object is masked from 'package:tidyr':
## 
##     smiths
## 
## Loading required package: ampvis2
## Loading required package: tsnemicrobiota

##       bookdown          knitr      tidyverse           plyr           grid 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##      gridExtra     kableExtra         xtable         ggpubr       phyloseq 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##          dada2       DECIPHER       phangorn         ggpubr    BiocManager 
##           TRUE           TRUE           TRUE           TRUE           TRUE 
##         DESeq2     microbiome          philr         btools      fantaxtic 
##           TRUE           TRUE          FALSE           TRUE           TRUE 
##        ampvis2 tsnemicrobiota 
##           TRUE           TRUE

#Extract sample names

miseq_path <- "./" # CHANGE to the directory containing the fastq files after unzipping.

# Sort ensures forward/reverse reads are in same order
fnFs <- sort(list.files(miseq_path, pattern="_R1_001.fastq"))

# Extract sample names, assuming filenames have format: SAMPLENAME_XXX.fastq
sampleNames <- sapply(strsplit(fnFs, "_"), `[`, 1)
# Specify the full path to the fnFs
fnFs <- file.path(miseq_path, fnFs)

fnFs[1:3]

## [1] ".//CUN02-V2_1_L001_R1_001.fastq" ".//CUN02-V3_4_L001_R1_001.fastq"
## [3] ".//CUN03-V2_7_L001_R1_001.fastq"

#plot quality profile two samples
plotQualityProfile(fnFs[1:2])

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

#plot quality profile 
plotQualityProfile(fnFs[c(1,11)])

## Warning: `guides(<scale> = FALSE)` is deprecated. Please use `guides(<scale> =
## "none")` instead.

#names(filtFs) <- sample.names
sample.names <- sapply(strsplit(basename(fnFs), "_"), `[`, 1)
head(sample.names)

## [1] "CUN02-V2" "CUN02-V3" "CUN03-V2" "CUN03-V3" "HCB01-V2" "HCB01-V3"

filt_path <- file.path(miseq_path, "filtered") # Place filtered files in filtered/ subdirectory
if(!file_test("-d", filt_path)) dir.create(filt_path)
filtFs <- file.path(filt_path, paste0(sampleNames, "_F_filt.fastq.gz"))

#Fiter and trim

outF1s<- filterAndTrim(fnFs, filtFs, truncLen=288, maxN=0, maxEE=2, truncQ=2, rm.phix=FALSE, compress=FALSE, verbose=TRUE)

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/CUN02-V2_F_filt.fastq.gz

## Read in 77582, output 55882 (72%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/CUN02-V3_F_filt.fastq.gz

## Read in 78500, output 56171 (71.6%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/CUN03-V2_F_filt.fastq.gz

## Read in 73044, output 50919 (69.7%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/CUN03-V3_F_filt.fastq.gz

## Read in 85148, output 60267 (70.8%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB01-V2_F_filt.fastq.gz

## Read in 77942, output 54420 (69.8%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB01-V3_F_filt.fastq.gz

## Read in 76092, output 51823 (68.1%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB02-V2_F_filt.fastq.gz

## Read in 62110, output 43706 (70.4%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB02-V3_F_filt.fastq.gz

## Read in 79907, output 56642 (70.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB03-V2_F_filt.fastq.gz

## Read in 90718, output 59814 (65.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB03-V3_F_filt.fastq.gz

## Read in 74075, output 49537 (66.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB04-V2_F_filt.fastq.gz

## Read in 78515, output 54934 (70%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB04-V3_F_filt.fastq.gz

## Read in 87173, output 60951 (69.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB06-V2_F_filt.fastq.gz

## Read in 86527, output 60214 (69.6%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB08-V2_F_filt.fastq.gz

## Read in 73143, output 49951 (68.3%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB09-V2_F_filt.fastq.gz

## Read in 92119, output 61560 (66.8%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB10-V2_F_filt.fastq.gz

## Read in 50022, output 35458 (70.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB12-V2_F_filt.fastq.gz

## Read in 73294, output 50478 (68.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HCB13-V3_F_filt.fastq.gz

## Read in 80869, output 55595 (68.7%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA01-V2_F_filt.fastq.gz

## Read in 67690, output 46666 (68.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA01-V3_F_filt.fastq.gz

## Read in 81275, output 57916 (71.3%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA02-V2_F_filt.fastq.gz

## Read in 81425, output 55645 (68.3%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA02-V3_F_filt.fastq.gz

## Read in 77772, output 53549 (68.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA03-V2_F_filt.fastq.gz

## Read in 78250, output 53339 (68.2%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA03-V3_F_filt.fastq.gz

## Read in 86178, output 60153 (69.8%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA04-V2_F_filt.fastq.gz

## Read in 103598, output 71280 (68.8%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA04-V3_F_filt.fastq.gz

## Read in 74184, output 50179 (67.6%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVA05-V3_F_filt.fastq.gz

## Read in 94023, output 63596 (67.6%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVR01-V2_F_filt.fastq.gz

## Read in 48149, output 32090 (66.6%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/HVR02-V2_F_filt.fastq.gz

## Read in 66957, output 44071 (65.8%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/ISBS01-V2_F_filt.fastq.gz

## Read in 82565, output 57700 (69.9%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/ISBS03-V2_F_filt.fastq.gz

## Read in 244723, output 168013 (68.7%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/ISBS06-V2_F_filt.fastq.gz

## Read in 81278, output 55876 (68.7%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/ISBS08-V2_F_filt.fastq.gz

## Read in 44945, output 28090 (62.5%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/ISBS09-V2_F_filt.fastq.gz

## Read in 60124, output 38568 (64.1%) filtered sequences.

## Overwriting file:/home/eduvb/qiime2/Rstudio/filtered/ISBS10-V2_F_filt.fastq.gz

## Read in 72508, output 48676 (67.1%) filtered sequences.

head(outF1s)

##                               reads.in reads.out
## CUN02-V2_1_L001_R1_001.fastq     77582     55882
## CUN02-V3_4_L001_R1_001.fastq     78500     56171
## CUN03-V2_7_L001_R1_001.fastq     73044     50919
## CUN03-V3_9_L001_R1_001.fastq     85148     60267
## HCB01-V2_13_L001_R1_001.fastq    77942     54420
## HCB01-V3_16_L001_R1_001.fastq    76092     51823

#Dereplication

derepFs <- derepFastq(filtFs, verbose=TRUE)

## Dereplicating sequence entries in Fastq file: .//filtered/CUN02-V2_F_filt.fastq.gz

## Encountered 19350 unique sequences from 55882 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/CUN02-V3_F_filt.fastq.gz

## Encountered 20061 unique sequences from 56171 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/CUN03-V2_F_filt.fastq.gz

## Encountered 20630 unique sequences from 50919 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/CUN03-V3_F_filt.fastq.gz

## Encountered 24207 unique sequences from 60267 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB01-V2_F_filt.fastq.gz

## Encountered 22104 unique sequences from 54420 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB01-V3_F_filt.fastq.gz

## Encountered 19675 unique sequences from 51823 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB02-V2_F_filt.fastq.gz

## Encountered 14022 unique sequences from 43706 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB02-V3_F_filt.fastq.gz

## Encountered 18830 unique sequences from 56642 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB03-V2_F_filt.fastq.gz

## Encountered 20981 unique sequences from 59814 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB03-V3_F_filt.fastq.gz

## Encountered 17912 unique sequences from 49537 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB04-V2_F_filt.fastq.gz

## Encountered 17659 unique sequences from 54934 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB04-V3_F_filt.fastq.gz

## Encountered 19467 unique sequences from 60951 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB06-V2_F_filt.fastq.gz

## Encountered 23786 unique sequences from 60214 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB08-V2_F_filt.fastq.gz

## Encountered 17432 unique sequences from 49951 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB09-V2_F_filt.fastq.gz

## Encountered 24622 unique sequences from 61560 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB10-V2_F_filt.fastq.gz

## Encountered 12603 unique sequences from 35458 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB12-V2_F_filt.fastq.gz

## Encountered 17518 unique sequences from 50478 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HCB13-V3_F_filt.fastq.gz

## Encountered 22507 unique sequences from 55595 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA01-V2_F_filt.fastq.gz

## Encountered 18268 unique sequences from 46666 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA01-V3_F_filt.fastq.gz

## Encountered 21055 unique sequences from 57916 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA02-V2_F_filt.fastq.gz

## Encountered 21588 unique sequences from 55645 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA02-V3_F_filt.fastq.gz

## Encountered 21872 unique sequences from 53549 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA03-V2_F_filt.fastq.gz

## Encountered 18099 unique sequences from 53339 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA03-V3_F_filt.fastq.gz

## Encountered 17148 unique sequences from 60153 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA04-V2_F_filt.fastq.gz

## Encountered 31215 unique sequences from 71280 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA04-V3_F_filt.fastq.gz

## Encountered 20455 unique sequences from 50179 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVA05-V3_F_filt.fastq.gz

## Encountered 25355 unique sequences from 63596 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVR01-V2_F_filt.fastq.gz

## Encountered 13935 unique sequences from 32090 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/HVR02-V2_F_filt.fastq.gz

## Encountered 17785 unique sequences from 44071 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/ISBS01-V2_F_filt.fastq.gz

## Encountered 21785 unique sequences from 57700 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/ISBS03-V2_F_filt.fastq.gz

## Encountered 56826 unique sequences from 168013 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/ISBS06-V2_F_filt.fastq.gz

## Encountered 25194 unique sequences from 55876 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/ISBS08-V2_F_filt.fastq.gz

## Encountered 12163 unique sequences from 28090 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/ISBS09-V2_F_filt.fastq.gz

## Encountered 13497 unique sequences from 38568 total sequences read.

## Dereplicating sequence entries in Fastq file: .//filtered/ISBS10-V2_F_filt.fastq.gz

## Encountered 17999 unique sequences from 48676 total sequences read.

# Name the derep-class objects by the sample names
names(derepFs) <- sampleNames

#Learn errors
errF <- learnErrors(filtFs, multithread=TRUE)

## 107478144 total bases in 373188 reads from 7 samples will be used for learning the error rates.

plotErrors(errF)

#DADA2

dadaFs <- dada(derepFs, err=errF, multithread=TRUE)

## Sample 1 - 55882 reads in 19350 unique sequences.
## Sample 2 - 56171 reads in 20061 unique sequences.
## Sample 3 - 50919 reads in 20630 unique sequences.
## Sample 4 - 60267 reads in 24207 unique sequences.
## Sample 5 - 54420 reads in 22104 unique sequences.
## Sample 6 - 51823 reads in 19675 unique sequences.
## Sample 7 - 43706 reads in 14022 unique sequences.
## Sample 8 - 56642 reads in 18830 unique sequences.
## Sample 9 - 59814 reads in 20981 unique sequences.
## Sample 10 - 49537 reads in 17912 unique sequences.
## Sample 11 - 54934 reads in 17659 unique sequences.
## Sample 12 - 60951 reads in 19467 unique sequences.
## Sample 13 - 60214 reads in 23786 unique sequences.
## Sample 14 - 49951 reads in 17432 unique sequences.
## Sample 15 - 61560 reads in 24622 unique sequences.
## Sample 16 - 35458 reads in 12603 unique sequences.
## Sample 17 - 50478 reads in 17518 unique sequences.
## Sample 18 - 55595 reads in 22507 unique sequences.
## Sample 19 - 46666 reads in 18268 unique sequences.
## Sample 20 - 57916 reads in 21055 unique sequences.
## Sample 21 - 55645 reads in 21588 unique sequences.
## Sample 22 - 53549 reads in 21872 unique sequences.
## Sample 23 - 53339 reads in 18099 unique sequences.
## Sample 24 - 60153 reads in 17148 unique sequences.
## Sample 25 - 71280 reads in 31215 unique sequences.
## Sample 26 - 50179 reads in 20455 unique sequences.
## Sample 27 - 63596 reads in 25355 unique sequences.
## Sample 28 - 32090 reads in 13935 unique sequences.
## Sample 29 - 44071 reads in 17785 unique sequences.
## Sample 30 - 57700 reads in 21785 unique sequences.
## Sample 31 - 168013 reads in 56826 unique sequences.
## Sample 32 - 55876 reads in 25194 unique sequences.
## Sample 33 - 28090 reads in 12163 unique sequences.
## Sample 34 - 38568 reads in 13497 unique sequences.
## Sample 35 - 48676 reads in 17999 unique sequences.

#Inspect dada2 object
dadaFs[[1]]

## dada-class: object describing DADA2 denoising results
## 383 sequence variants were inferred from 19350 input unique sequences.
## Key parameters: OMEGA_A = 1e-40, OMEGA_C = 1e-40, BAND_SIZE = 16

##Construct sequence table and remove chimeras

library(dada2)
library(Rcpp)
#BiocManager::install('taRifx')
library(taRifx)

## 
## Attaching package: 'taRifx'

## The following object is masked from 'package:SummarizedExperiment':
## 
##     shift

## The following object is masked from 'package:GenomicRanges':
## 
##     shift

## The following object is masked from 'package:IRanges':
## 
##     shift

## The following object is masked from 'package:S4Vectors':
## 
##     first

## The following objects are masked from 'package:dplyr':
## 
##     between, distinct, first, last

## The following object is masked from 'package:purrr':
## 
##     rep_along

mergers <-merge.list(dadaFs, derepFs)

seqtabAll <- makeSequenceTable(mergers[!grepl("Mock", names(mergers))])
table(nchar(getSequences(seqtabAll)))

## 
##   288 
## 10484

#remove chimeric sequences 
seqtabNoC <- removeBimeraDenovo(seqtabAll)

#Assign taxonomy <<>>

#SILVA Reference Genetic Database

#wget \
#  -O "silva_nr99_v138.1_train_set.fa.gz" \
#  "https://zenodo.org/record/4587955/files/silva_nr99_v138.1_train_set.fa.gz?download=1"

#wget \
#  -O "silva_species_assignment_v138.1.fa.gz" \
#  "https://zenodo.org/record/4587955/files/silva_species_assignment_v138.1.fa.gz?download=1"

tt <- assignTaxonomy(seqtabNoC, "silva_nr99_v138.1_train_set.fa.gz")
#tt.plus <- addSpecies(tt, "silva_species_assignment_v138.1.fa.gz", verbose=TRUE)
#load("tt.plus.Rda")

#unname(head(tt.plus))
unname(head(tt))

##      [,1]       [,2]             [,3]                  [,4]              
## [1,] "Bacteria" "Firmicutes"     "Bacilli"             "Lactobacillales" 
## [2,] "Bacteria" "Firmicutes"     "Bacilli"             "Lactobacillales" 
## [3,] "Bacteria" "Firmicutes"     "Clostridia"          "Lachnospirales"  
## [4,] "Bacteria" "Firmicutes"     "Bacilli"             "Lactobacillales" 
## [5,] "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Enterobacterales"
## [6,] "Bacteria" "Proteobacteria" "Gammaproteobacteria" "Enterobacterales"
##      [,5]                 [,6]                  
## [1,] "Streptococcaceae"   "Streptococcus"       
## [2,] "Streptococcaceae"   "Streptococcus"       
## [3,] "Lachnospiraceae"    "Blautia"             
## [4,] "Streptococcaceae"   "Streptococcus"       
## [5,] "Enterobacteriaceae" "Escherichia-Shigella"
## [6,] "Enterobacteriaceae" "Escherichia-Shigella"

taxTab <- tt

#Construct phylogenetic tree

#using the DECIPHER R package
seqs <- getSequences(seqtabNoC)
names(seqs) <- seqs # This propagates to the tip labels of the tree
alignment <- AlignSeqs(DNAStringSet(seqs), anchor=NA,verbose=FALSE)


#phangorn R package is then used to construct a phylogenetic tree

phangAlign <- phyDat(as(alignment, "matrix"), type="DNA")
dm <- dist.ml(phangAlign)

treeNJ <- NJ(dm) # Note, tip order != sequence order

fit = pml(treeNJ, data=phangAlign)

## negative edges length changed to 0!

fitGTR <- update(fit, k=4, inv=0.2)
fitGTR <- optim.pml(fitGTR, model="GTR", optInv=TRUE, optGamma=TRUE,
        rearrangement = "stochastic", control = pml.control(trace = 0))
detach("package:phangorn", unload=TRUE)

##Combine data into a phyloseq object

#Load metadata file
samdf <- read.csv("Metadata_OTUs_mod_1.csv",header=TRUE)
rownames(samdf) <- samdf$SampleID

keep.cols <- c("SampleID",  "Edad_a_la_inclusion", "Sexo",  "Respuesta_100",    "aloTPH_previo_si_no",  "Citogenética_alto_riesgo", "Afect_Extramedular_PET_screening_si_no",   "Dias_producción",  "CRS_si_no", "Grado_máx_CRS",   "HOSPITAL", "Visit")

samdf <- samdf[rownames(seqtabAll), keep.cols]

# the sample metadata, the sequence taxonomies, and the phylogenetic tree – can now be combined into a single object
ps <- phyloseq(otu_table(seqtabNoC, taxa_are_rows=FALSE), 
               sample_data(samdf), 
               tax_table(taxTab),phy_tree(fitGTR$tree))

ps

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 6149 taxa and 35 samples ]
## sample_data() Sample Data:       [ 35 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 6149 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 6149 tips and 6147 internal nodes ]

#Loading required libraries

library(tidyverse)
#BiocManager::install("curatedMetagenomicData")
library(curatedMetagenomicData)

## Loading required package: AnnotationHub

## Loading required package: BiocFileCache

## Loading required package: dbplyr

## 
## Attaching package: 'dbplyr'

## The following objects are masked from 'package:dplyr':
## 
##     ident, sql

## 
## Attaching package: 'AnnotationHub'

## The following object is masked from 'package:Biobase':
## 
##     cache

## Loading required package: ExperimentHub

##                       _           _
##   ___ _   _ _ __ __ _| |_ ___  __| |
##  / __| | | | '__/ _` | __/ _ \/ _` |
## | (__| |_| | | | (_| | ||  __/ (_| |
##  \___|\__,_|_|  \__,_|\__\___|\__,_|
##  __  __      _                                         _
## |  \/  | ___| |_ __ _  __ _  ___ _ __   ___  _ __ ___ (_) ___
## | |\/| |/ _ \ __/ _` |/ _` |/ _ \ '_ \ / _ \| '_ ` _ \| |/ __|
## | |  | |  __/ || (_| | (_| |  __/ | | | (_) | | | | | | | (__
## |_|  |_|\___|\__\__,_|\__, |\___|_| |_|\___/|_| |_| |_|_|\___|
##  ____        _        |___/
## |  _ \  __ _| |_ __ _
## | | | |/ _` | __/ _` |
## | |_| | (_| | || (_| |
## |____/ \__,_|\__\__,_|

library(phyloseq)
library(DESeq2)
#BiocManager::install("apeglm")
library(apeglm)
#BiocManager::install("scran")
#BiocManager::install("zinbwave")
library(zinbwave)

## Loading required package: SingleCellExperiment

library(scran)
library(cowplot)

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggpubr':
## 
##     get_legend

#BiocManager::install("VennDiagram")
library(VennDiagram)

## Loading required package: futile.logger

## 
## Attaching package: 'VennDiagram'

## The following object is masked from 'package:ape':
## 
##     rotate

## The following object is masked from 'package:ggpubr':
## 
##     rotate

#BiocManager::install("pscl")
library(pscl)

## Classes and Methods for R developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University
## Simon Jackman
## hurdle and zeroinfl functions by Achim Zeileis

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

#####Phyloseq object#####

#Taxonomic Filtering
# Show available ranks in the dataset
rank_names(ps)

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"

## Create table, number of features for each phyla
table(tax_table(ps)[, "Phylum"], exclude = NULL)

## 
##  Actinobacteriota      Bacteroidota  Campylobacterota     Cyanobacteria 
##               622               749                 6                 9 
##  Desulfobacterota     Euryarchaeota        Firmicutes    Fusobacteriota 
##                21                 8              4290                12 
##   Patescibacteria    Proteobacteria      Synergistota Verrucomicrobiota 
##                18               246                12                30 
##              <NA> 
##               126

#features were annotated with a Phylum of NA. These features are probably artifacts in a dataset like this, and should be removed.
ps <- subset_taxa(ps, !is.na(Phylum) & !Phylum %in% c("", "uncharacterized"))

# Compute prevalence of each feature, store as data.frame
prevdf = apply(X = otu_table(ps),
               MARGIN = ifelse(taxa_are_rows(ps), yes = 1, no = 2),
               FUN = function(x){sum(x > 0)})
# Add taxonomy and total read counts to this data.frame
prevdf = data.frame(Prevalence = prevdf,
                    TotalAbundance = taxa_sums(ps),
                    tax_table(ps))

#Create a new phyloseq-object
psv2 <- ps
#Define levels as factor
sample_data(psv2)$Respuesta_100 <- as.factor(sample_data(psv2)$Respuesta_100)
levels(sample_data(psv2)$Respuesta_100)

## [1] "PR"   "sCR"  "VGPR"

#rename <<VGPR>> en <<sCR>>
levels(sample_data(psv2)$Respuesta_100)[3] <- "sCR"
levels(sample_data(psv2)$Respuesta_100)

## [1] "PR"  "sCR"

sample_data(psv2)$Respuesta_100

##  [1] sCR sCR sCR sCR PR  PR  PR  PR  sCR sCR sCR sCR sCR PR  sCR sCR sCR sCR PR 
## [20] PR  sCR sCR sCR sCR sCR sCR PR  sCR sCR sCR sCR PR  sCR PR  sCR
## Levels: PR sCR

#Create a new phyloseq-object filtered by <<Visit==V2>>

psv2 <- subset_samples(psv2, Visit=="V2" )
sample_data(psv2)

##            SampleID Edad_a_la_inclusion   Sexo Respuesta_100
## CUN02-V2   CUN02-V2                  43  MUJER           sCR
## CUN03-V2   CUN03-V2                  45 HOMBRE           sCR
## HCB01-V2   HCB01-V2                  66 HOMBRE            PR
## HCB02-V2   HCB02-V2                  66 HOMBRE            PR
## HCB03-V2   HCB03-V2                  60 HOMBRE           sCR
## HCB04-V2   HCB04-V2                  55  MUJER           sCR
## HCB06-V2   HCB06-V2                  73  MUJER           sCR
## HCB08-V2   HCB08-V2                  62  MUJER            PR
## HCB09-V2   HCB09-V2                  53 HOMBRE           sCR
## HCB10-V2   HCB10-V2                  52 HOMBRE           sCR
## HCB12-V2   HCB12-V2                  64 HOMBRE           sCR
## HVA01-V2   HVA01-V2                  58 HOMBRE            PR
## HVA02-V2   HVA02-V2                  67 HOMBRE           sCR
## HVA03-V2   HVA03-V2                  61 HOMBRE           sCR
## HVA04-V2   HVA04-V2                  74  MUJER           sCR
## HVR01-V2   HVR01-V2                  65 HOMBRE           sCR
## HVR02-V2   HVR02-V2                  46 HOMBRE           sCR
## ISBS01-V2 ISBS01-V2                  56 HOMBRE           sCR
## ISBS03-V2 ISBS03-V2                  70  MUJER           sCR
## ISBS06-V2 ISBS06-V2                  60 HOMBRE            PR
## ISBS08-V2 ISBS08-V2                  53 HOMBRE           sCR
## ISBS09-V2 ISBS09-V2                  65  MUJER            PR
## ISBS10-V2 ISBS10-V2                  68 HOMBRE           sCR
##           aloTPH_previo_si_no Citogenética_alto_riesgo
## CUN02-V2                   NO                       NO
## CUN03-V2                   NO                       SI
## HCB01-V2                   NO                       NO
## HCB02-V2                   NO                       SI
## HCB03-V2                   NO                       SI
## HCB04-V2                   NO                       NO
## HCB06-V2                   NO                       NO
## HCB08-V2                   NO                       NO
## HCB09-V2                   NO                       NO
## HCB10-V2                   SI                       NO
## HCB12-V2                   NO                       SI
## HVA01-V2                   NO                       NO
## HVA02-V2                   NO                       NO
## HVA03-V2                   SI                       NO
## HVA04-V2                   NO                       NO
## HVR01-V2                                            NO
## HVR02-V2                                            NO
## ISBS01-V2                  NO                       SI
## ISBS03-V2                  NO                         
## ISBS06-V2                  NO                       NO
## ISBS08-V2                  NO                         
## ISBS09-V2                  NO                       NO
## ISBS10-V2                  NO                       SI
##           Afect_Extramedular_PET_screening_si_no Dias_producción CRS_si_no
## CUN02-V2                                      SI             120        SI
## CUN03-V2                                      NO             100        SI
## HCB01-V2                                      SI             110        SI
## HCB02-V2                                      NO             110        SI
## HCB03-V2                                      SI             110        SI
## HCB04-V2                                      NO             110        SI
## HCB06-V2                                      NO              90        SI
## HCB08-V2                                      SI             100        SI
## HCB09-V2                                      NO              90        SI
## HCB10-V2                                      SI             100        SI
## HCB12-V2                                      NO             110        NO
## HVA01-V2                                      NO             120        SI
## HVA02-V2                                      SI             100        SI
## HVA03-V2                                      SI             140        SI
## HVA04-V2                                      NO             100        SI
## HVR01-V2                                      SI             110        SI
## HVR02-V2                                      SI             120        SI
## ISBS01-V2                                     SI             120        SI
## ISBS03-V2                                     SI             100        SI
## ISBS06-V2                                     NO             130        SI
## ISBS08-V2                                     SI             110        SI
## ISBS09-V2                                     SI             100        SI
## ISBS10-V2                                     SI             140        NO
##           Grado_máx_CRS HOSPITAL Visit
## CUN02-V2             10      CUN    V2
## CUN03-V2             20      CUN    V2
## HCB01-V2             20      HCB    V2
## HCB02-V2             20      HCB    V2
## HCB03-V2             10      HCB    V2
## HCB04-V2             10      HCB    V2
## HCB06-V2             10      HCB    V2
## HCB08-V2             10      HCB    V2
## HCB09-V2             10      HCB    V2
## HCB10-V2             10      HCB    V2
## HCB12-V2             NA      HCB    V2
## HVA01-V2             20      HVA    V2
## HVA02-V2             20      HVA    V2
## HVA03-V2             10      HVA    V2
## HVA04-V2             10      HVA    V2
## HVR01-V2             10      HVR    V2
## HVR02-V2             10      HVR    V2
## ISBS01-V2            10     ISBS    V2
## ISBS03-V2            10     ISBS    V2
## ISBS06-V2            10     ISBS    V2
## ISBS08-V2            20     ISBS    V2
## ISBS09-V2            10     ISBS    V2
## ISBS10-V2            NA     ISBS    V2

#Review the levels of the variable  <<Respuesta_100>>
levels(sample_data(psv2)$Respuesta_100)

## [1] "PR"  "sCR"

#Summary variable <<Respuesta_100>>
summary(sample_data(psv2)$Respuesta_100)

##  PR sCR 
##   6  17

list(sample_data(psv2)$Visit)

## [[1]]
##  [1] "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2"
## [16] "V2" "V2" "V2" "V2" "V2" "V2" "V2" "V2"

#Filtered phyloseq object by patients visit <<V2>>
psv2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 6023 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 6023 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 6023 tips and 6021 internal nodes ]

#Taxonomic prevalence

# Show available ranks in the dataset
rank_names(psv2)

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"

## Create table, number of features for each phyla
table(tax_table(psv2)[, "Phylum"], exclude = NULL)

## 
##  Actinobacteriota      Bacteroidota  Campylobacterota     Cyanobacteria 
##               622               749                 6                 9 
##  Desulfobacterota     Euryarchaeota        Firmicutes    Fusobacteriota 
##                21                 8              4290                12 
##   Patescibacteria    Proteobacteria      Synergistota Verrucomicrobiota 
##                18               246                12                30

#features were annotated with a Phylum of NA. These features are probably artifacts in a dataset like this, and should be removed.
psv2 <- subset_taxa(psv2, !is.na(Phylum) & !Phylum %in% c("", "uncharacterized"))


# compute prevalence for each feature and store it in a data frame
prevdf = apply(X = otu_table(psv2),
               MARGIN = ifelse(taxa_are_rows(psv2), yes = 1, no = 2),
               FUN = function(x){sum(x > 0)})
# add the taxonomy
prevdf = data.frame(Prevalence = prevdf,
                    TotalAbundance = taxa_sums(psv2),
                    tax_table(psv2))

plyr::ddply(prevdf, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))}) -> dfprev
kable(dfprev)

Phylum	1	2
Actinobacteriota	1.1752412	731
Bacteroidota	1.2616822	945
Campylobacterota	0.8333333	5
Cyanobacteria	0.8888889	8
Desulfobacterota	1.1904762	25
Euryarchaeota	2.5000000	20
Firmicutes	1.2559441	5388
Fusobacteriota	0.5000000	6
Patescibacteria	0.8333333	15
Proteobacteria	1.1260163	277
Synergistota	1.0000000	12
Verrucomicrobiota	2.2666667	68

# Show available ranks in the dataset
rank_names(psv2)

## [1] "Kingdom" "Phylum"  "Class"   "Order"   "Family"  "Genus"

## Create table, number of features for each phyla
table(tax_table(psv2)[, "Phylum"], exclude = NULL)

## 
##  Actinobacteriota      Bacteroidota  Campylobacterota     Cyanobacteria 
##               622               749                 6                 9 
##  Desulfobacterota     Euryarchaeota        Firmicutes    Fusobacteriota 
##                21                 8              4290                12 
##   Patescibacteria    Proteobacteria      Synergistota Verrucomicrobiota 
##                18               246                12                30

#features were annotated with a Phylum of NA. These features are probably artifacts in a dataset like this, and should be removed.
psv2 <- subset_taxa(psv2, !is.na(Phylum) & !Phylum %in% c("", "uncharacterized"))


# Computamos prevalencia para cada feature y la guardamos en un data frame
prevdf = apply(X = otu_table(psv2),
               MARGIN = ifelse(taxa_are_rows(psv2), yes = 1, no = 2),
               FUN = function(x){sum(x > 0)})
# Le agregamos la taxonomía
prevdf = data.frame(Prevalence = prevdf,
                    TotalAbundance = taxa_sums(psv2),
                    tax_table(psv2))

plyr::ddply(prevdf, "Phylum", function(df1){cbind(mean(df1$Prevalence),sum(df1$Prevalence))}) -> dfprev
kable(dfprev)

Phylum	1	2
Actinobacteriota	1.1752412	731
Bacteroidota	1.2616822	945
Campylobacterota	0.8333333	5
Cyanobacteria	0.8888889	8
Desulfobacterota	1.1904762	25
Euryarchaeota	2.5000000	20
Firmicutes	1.2559441	5388
Fusobacteriota	0.5000000	6
Patescibacteria	0.8333333	15
Proteobacteria	1.1260163	277
Synergistota	1.0000000	12
Verrucomicrobiota	2.2666667	68

#La columna 1 representa la media de read counts para ese Phylum, mientras que la columna 2 representa la suma

# Define phyla to filter
filterPhyla = c("Campylobacterota", "Cyanobacteria", "Halanaerobiaeota", "Fusobacteriota", "Synergistota" )
psv2 = subset_taxa(psv2, !Phylum %in% filterPhyla)
psv2

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 5984 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 5984 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 5984 tips and 5982 internal nodes ]

#Filter taxa

# ilter taxa according to a threshold of mean number of _read counts_, in this case 1e-5
psd2 <- filter_taxa(psv2, function(x) mean(x) > 1e-5, TRUE)

# remove taxa that are not observed more than X times in at least 10% of the samples
psd3 <- filter_taxa(psd2, function(x) sum(x > 2) > (0.1*length(x)), TRUE)

# filter samples with less than 1000 reads
psd4 = prune_samples(sample_sums(psd3) > 1000, psd3)

psd4

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 557 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 557 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 557 tips and 555 internal nodes ]

# We select the taxa of interest (ALLS)
prevdf1 = subset(prevdf, Phylum %in% get_taxa_unique(psd4, "Phylum"))
ggplot(prevdf1, aes(TotalAbundance, Prevalence / nsamples(psv2),color=Phylum)) +
  # Agregamos una línea para nuestro umbral
  geom_hline(yintercept = 0.05, alpha = 0.5, linetype = 2) +  geom_point(size = 2, alpha = 0.7) +
  scale_x_log10() +  xlab("Total Abundance") + ylab("Prevalence [Frac. Samples]") +
  facet_wrap(~Phylum) + theme(legend.position="none")

## Warning: Transformation introduced infinite values in continuous x-axis

# replace the sequences with a generic name
psd5 <- psd4
taxa_names(psd5) <- paste0("ASV", seq(ntaxa(psd5)))
head(taxa_names(psd5))

## [1] "ASV1" "ASV2" "ASV3" "ASV4" "ASV5" "ASV6"

#plot the distribution of read counts by sample

sample_sum_df <- data.frame(sum = sample_sums(psd5))

ggplot(sample_sum_df, aes(x = sum)) + 
  geom_histogram(color = "black", fill = "grey", binwidth = 2500) +
  ggtitle("Distribution of sample sequencing depth") + 
  xlab("Read counts") +
  theme(axis.title.y = element_blank())

#calculate rarefaction curves for each sample

# First  load  scripts 
scripts <- c("graphical_methods.R",
             "tree_methods.R",
             "plot_merged_trees.R",
             "specificity_methods.R",
             "ternary_plot.R",
             "richness.R",
             "edgePCA.R",
             "copy_number_correction.R",
             "import_frogs.R",
             "prevalence.R",
             "compute_niche.R")
urls <- paste0("https://raw.githubusercontent.com/mahendra-mariadassou/phyloseq-extended/master/R/", scripts)

for (url in urls) {
  source(url)
}

## Loading required package: scales

## 
## Attaching package: 'scales'

## The following object is masked from 'package:microbiome':
## 
##     alpha

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

## Loading required package: vegan

## Loading required package: permute

## Warning: package 'permute' was built under R version 4.0.5

## Loading required package: lattice

## This is vegan 2.5-7

## 
## Attaching package: 'vegan'

## The following object is masked from 'package:microbiome':
## 
##     diversity

#Plot the rarefaction curves

p <- ggrare(psd5, step = 100, color = "Respuesta_100", label = "SampleID", se = TRUE)

## rarefying sample CUN02-V2
## rarefying sample CUN03-V2
## rarefying sample HCB01-V2
## rarefying sample HCB02-V2
## rarefying sample HCB03-V2
## rarefying sample HCB04-V2
## rarefying sample HCB06-V2
## rarefying sample HCB08-V2
## rarefying sample HCB09-V2
## rarefying sample HCB10-V2
## rarefying sample HCB12-V2
## rarefying sample HVA01-V2
## rarefying sample HVA02-V2
## rarefying sample HVA03-V2
## rarefying sample HVA04-V2
## rarefying sample HVR01-V2
## rarefying sample HVR02-V2
## rarefying sample ISBS01-V2
## rarefying sample ISBS03-V2
## rarefying sample ISBS06-V2
## rarefying sample ISBS08-V2
## rarefying sample ISBS09-V2
## rarefying sample ISBS10-V2

(p <- p + facet_wrap(~Respuesta_100))

Table of read counts

#For visualization in HTML format
library("kableExtra")

head(otu_table(psd5)) %>%
  kable(format = "html", col.names = colnames(otu_table(psd5))) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "350px")

	ASV1	ASV2	ASV3	ASV4	ASV5	ASV6	ASV7	ASV8	ASV9	ASV10	ASV11	ASV12	ASV13	ASV14	ASV15	ASV16	ASV17	ASV18	ASV19	ASV20	ASV21	ASV22	ASV23	ASV24	ASV25	ASV26	ASV27	ASV28	ASV30	ASV31	ASV32	ASV33	ASV34	ASV35	ASV36	ASV37	ASV38	ASV39	ASV41	ASV42	ASV43	ASV44	ASV45	ASV47	ASV48	ASV49	ASV50	ASV51	ASV52	ASV53	ASV54	ASV55	ASV56	ASV57	ASV58	ASV59	ASV60	ASV61	ASV62	ASV63	ASV64	ASV65	ASV66	ASV67	ASV68	ASV69	ASV70	ASV71	ASV72	ASV73	ASV74	ASV75	ASV76	ASV77	ASV78	ASV79	ASV80	ASV81	ASV82	ASV83	ASV84	ASV85	ASV86	ASV87	ASV88	ASV89	ASV90	ASV91	ASV92	ASV93	ASV94	ASV95	ASV96	ASV97	ASV99	ASV100	ASV101	ASV102	ASV103	ASV104	ASV105	ASV106	ASV107	ASV108	ASV109	ASV110	ASV111	ASV112	ASV113	ASV114	ASV115	ASV116	ASV117	ASV119	ASV120	ASV121	ASV122	ASV123	ASV124	ASV125	ASV126	ASV127	ASV128	ASV129	ASV130	ASV131	ASV132	ASV133	ASV134	ASV135	ASV136	ASV137	ASV138	ASV139	ASV140	ASV141	ASV142	ASV143	ASV144	ASV145	ASV147	ASV148	ASV149	ASV151	ASV153	ASV154	ASV155	ASV156	ASV157	ASV158	ASV159	ASV160	ASV161	ASV163	ASV164	ASV165	ASV166	ASV167	ASV169	ASV170	ASV171	ASV172	ASV173	ASV175	ASV176	ASV177	ASV178	ASV179	ASV180	ASV183	ASV184	ASV185	ASV186	ASV187	ASV188	ASV189	ASV190	ASV191	ASV192	ASV193	ASV194	ASV195	ASV196	ASV197	ASV198	ASV200	ASV201	ASV202	ASV203	ASV204	ASV205	ASV206	ASV207	ASV209	ASV210	ASV211	ASV212	ASV213	ASV214	ASV215	ASV216	ASV217	ASV221	ASV222	ASV223	ASV224	ASV225	ASV227	ASV228	ASV230	ASV233	ASV234	ASV235	ASV236	ASV237	ASV238	ASV239	ASV240	ASV241	ASV242	ASV243	ASV244	ASV246	ASV248	ASV249	ASV250	ASV251	ASV252	ASV253	ASV254	ASV255	ASV256	ASV257	ASV258	ASV259	ASV260	ASV261	ASV263	ASV264	ASV266	ASV267	ASV268	ASV269	ASV270	ASV271	ASV273	ASV274	ASV275	ASV277	ASV278	ASV279	ASV281	ASV282	ASV283	ASV284	ASV285	ASV287	ASV288	ASV289	ASV290	ASV291	ASV292	ASV293	ASV294	ASV295	ASV296	ASV300	ASV302	ASV303	ASV304	ASV305	ASV306	ASV307	ASV308	ASV309	ASV310	ASV311	ASV312	ASV313	ASV314	ASV315	ASV316	ASV317	ASV318	ASV319	ASV321	ASV322	ASV323	ASV324	ASV325	ASV326	ASV327	ASV328	ASV329	ASV330	ASV331	ASV332	ASV333	ASV335	ASV337	ASV338	ASV339	ASV340	ASV341	ASV342	ASV343	ASV344	ASV346	ASV347	ASV348	ASV349	ASV350	ASV351	ASV352	ASV353	ASV355	ASV356	ASV357	ASV359	ASV360	ASV361	ASV362	ASV363	ASV366	ASV367	ASV368	ASV370	ASV371	ASV372	ASV375	ASV376	ASV377	ASV378	ASV380	ASV383	ASV384	ASV385	ASV386	ASV387	ASV389	ASV390	ASV392	ASV394	ASV395	ASV396	ASV399	ASV400	ASV401	ASV403	ASV404	ASV405	ASV407	ASV409	ASV410	ASV416	ASV418	ASV422	ASV423	ASV424	ASV426	ASV429	ASV430	ASV432	ASV434	ASV435	ASV437	ASV438	ASV439	ASV440	ASV444	ASV453	ASV455	ASV457	ASV459	ASV460	ASV461	ASV467	ASV470	ASV471	ASV472	ASV474	ASV475	ASV477	ASV479	ASV480	ASV481	ASV482	ASV483	ASV486	ASV487	ASV490	ASV492	ASV493	ASV495	ASV496	ASV497	ASV500	ASV501	ASV503	ASV505	ASV507	ASV508	ASV510	ASV511	ASV512	ASV514	ASV515	ASV516	ASV518	ASV519	ASV520	ASV522	ASV523	ASV526	ASV529	ASV531	ASV532	ASV534	ASV535	ASV537	ASV538	ASV541	ASV545	ASV546	ASV547	ASV549	ASV552	ASV553	ASV554	ASV555	ASV556	ASV557
CUN02-V2	0	0	528	0	0	0	447	1933	0	0	1863	272	0	213	0	0	1121	0	0	0	747	0	0	0	0	178	0	0	0	0	0	0	0	0	0	0	617	0	0	0	771	297	547	0	220	0	158	0	0	0	0	0	0	1099	718	334	110	188	0	0	0	0	0	745	0	0	0	0	0	0	0	0	0	0	0	380	0	0	0	0	0	0	0	1020	397	0	120	0	264	0	128	135	0	0	0	0	1983	2081	0	127	0	449	0	0	0	0	0	490	0	0	0	0	0	0	121	445	106	643	0	0	0	0	85	0	484	0	0	0	375	12	0	0	0	0	0	0	0	0	303	0	117	79	381	0	0	0	42	75	0	143	0	459	1247	0	357	388	68	0	0	0	225	0	395	0	0	0	0	0	245	135	0	0	0	341	79	0	0	0	57	0	0	0	0	0	221	0	185	852	215	0	93	0	157	0	0	0	0	0	0	59	0	0	0	52	0	0	0	0	204	90	44	0	118	217	0	0	0	0	0	104	0	0	28	0	0	0	0	0	0	0	0	88	0	0	91	0	0	0	0	0	283	0	0	0	0	0	0	0	56	0	0	0	0	125	108	0	0	104	0	269	0	0	0	0	0	0	0	0	84	0	0	267	0	69	0	49	591	111	0	41	0	0	0	0	0	89	0	0	0	0	134	0	53	0	0	0	0	0	51	141	0	0	0	0	0	0	0	0	0	0	0	0	274	104	125	0	0	0	0	181	0	0	0	0	0	44	176	171	125	47	0	0	112	0	0	0	0	0	0	103	0	0	0	0	0	0	0	42	0	60	0	0	0	0	0	144	0	0	149	0	0	0	87	0	108	0	0	0	0	52	0	6	0	0	61	168	0	0	0	0	0	0	0	0	0	9	0	0	0	0	0	0	0	0	0	0	0	0	0	59	0	0	0	46	0	87	0	0	0	59	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0
CUN03-V2	243	241	183	183	481	454	159	0	141	308	0	142	737	86	230	170	0	110	710	133	0	100	492	104	125	170	189	1340	312	432	83	76	1391	111	69	327	78	23	0	90	0	141	44	0	0	65	0	205	130	924	0	0	205	94	141	0	0	0	0	0	0	95	212	135	13	0	0	0	158	0	184	152	0	0	0	0	764	60	186	0	197	0	184	97	0	0	0	0	29	0	110	0	117	64	0	0	0	0	166	51	0	91	0	0	0	0	0	110	139	100	0	0	69	60	0	0	0	82	50	114	569	0	0	182	130	0	150	91	86	0	124	33	0	0	0	60	107	0	91	57	137	0	0	73	168	0	33	0	88	42	0	50	0	56	186	227	0	0	0	0	177	0	269	0	0	146	156	0	0	87	130	0	0	207	0	51	0	0	0	0	106	76	0	370	216	0	0	0	113	302	91	0	0	12	0	0	0	0	0	108	0	60	0	0	0	161	0	154	197	0	0	0	0	134	61	0	0	0	0	0	0	113	90	0	123	0	0	0	0	0	0	96	139	0	0	0	157	0	0	0	0	0	202	0	0	0	0	0	59	0	0	0	0	104	0	58	46	108	0	0	83	0	160	0	123	0	0	0	0	102	0	0	0	0	38	0	0	0	0	0	0	0	0	97	13	0	113	68	86	0	0	0	74	0	0	0	0	145	106	0	0	0	115	0	57	0	0	0	83	58	0	83	0	0	0	85	0	12	84	0	0	0	0	0	0	0	0	0	0	0	61	0	0	180	59	0	0	0	57	0	0	178	146	0	0	0	0	0	187	141	0	83	27	0	0	0	52	0	0	32	37	149	0	97	0	93	0	114	0	0	0	0	119	51	0	0	0	0	0	61	0	94	0	0	0	0	33	123	52	0	39	0	0	0	0	0	0	71	0	0	83	152	0	0	158	0	0	7	0	0	0	0	17	54	3	51	60	0	39	0	0	0	0	0	0	0	34	0	7	0	0	0
HCB01-V2	131	136	408	119	0	0	372	0	59	0	0	271	0	220	0	300	0	58	0	352	0	363	0	215	310	247	0	1192	1145	0	170	0	1185	198	0	1078	0	0	0	535	0	223	0	0	0	193	182	774	498	854	0	0	229	279	107	0	88	0	0	638	249	0	248	90	617	187	0	55	83	79	221	205	0	0	0	0	706	406	154	0	634	0	218	225	0	0	214	0	0	0	117	0	139	296	223	0	0	0	114	220	0	66	0	0	0	0	0	0	87	0	0	376	141	0	144	89	168	179	0	0	449	0	0	114	0	190	149	140	63	122	444	110	0	0	0	0	81	0	0	0	0	139	0	106	0	0	62	71	82	175	155	139	0	0	0	0	115	89	141	130	0	0	0	0	0	107	78	0	0	0	0	0	90	0	0	0	79	0	114	0	0	0	101	303	0	0	82	0	0	310	0	316	172	0	0	64	0	0	304	0	0	41	294	221	0	0	101	0	0	223	0	0	0	0	0	76	77	0	221	0	284	0	0	0	0	337	224	140	0	38	313	65	0	96	0	0	0	204	60	29	0	0	217	170	0	142	0	115	0	97	0	153	68	0	48	0	218	0	97	0	0	159	0	117	0	166	107	0	0	273	137	0	63	0	0	0	0	0	142	136	218	196	98	0	160	0	0	0	108	86	0	0	0	42	0	33	141	0	0	0	151	74	0	54	0	0	0	0	0	0	60	0	0	0	0	0	267	0	0	0	73	95	27	0	36	0	0	0	0	0	0	102	0	0	0	117	78	0	0	127	0	0	87	87	227	0	43	0	0	0	0	109	97	0	57	0	0	0	0	100	0	0	0	0	0	0	62	0	30	0	50	0	85	0	0	0	35	0	71	0	0	0	0	95	0	0	0	0	0	0	0	102	31	110	0	0	43	0	0	0	0	41	0	0	0	0	0	0	61	35	35	71	0	0	67	0	75	81	37	25	10	0	0	15	13	0	0	0	0	6	0	6
HCB02-V2	250	0	0	0	0	0	0	413	54	0	396	62	0	0	0	0	277	0	0	0	244	0	0	0	0	0	0	0	124	0	0	0	0	39	0	0	0	0	0	0	177	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	273	308	548	114	0	0	0	0	0	112	0	0	563	0	28	0	0	0	106	86	191	0	0	0	0	0	0	0	0	0	0	0	0	196	0	0	0	0	207	0	0	0	349	0	0	0	0	0	0	62	100	0	0	0	0	0	0	0	0	0	0	0	0	135	0	0	127	0	0	0	259	0	0	0	0	101	64	0	0	0	0	239	0	0	0	0	0	0	457	0	0	0	80	0	0	249	0	0	0	0	0	1639	0	0	0	0	0	427	0	0	0	0	0	1462	120	0	0	0	0	0	51	0	0	0	54	0	116	0	15	0	0	179	0	0	105	0	0	152	0	0	0	0	0	81	0	0	90	304	116	0	0	0	0	0	0	0	112	879	0	29	0	0	138	0	0	0	0	129	0	0	0	88	0	0	0	0	0	0	0	196	0	0	0	0	0	0	0	0	0	0	0	0	0	0	58	0	97	0	0	0	0	65	0	0	0	0	0	62	0	93	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	82	0	170	0	0	0	0	58	0	0	0	84	62	77	0	117	138	0	0	0	0	0	0	0	0	0	0	0	0	32	0	34	0	0	0	0	0	0	0	0	0	0	0	0	0	60	96	0	0	0	0	0	0	0	129	0	57	0	24	0	90	0	0	0	0	0	97	124	0	0	0	0	0	0	0	0	0	0	0	0	86	0	0	0	0	0	0	22	0	0	0	0	0	0	0	0	0	57	0	0	0	52	151	0	0	0	0	0	8	0	0	0	75	0	0	0	0	0	55	0	0	0	0	0	0	0	0	0	24	0	0	0	0	0	0	0	0	0	0
HCB03-V2	0	0	664	0	92	50	559	0	0	63	0	678	385	481	53	165	0	0	326	138	0	73	320	170	0	0	0	0	364	305	187	0	0	0	0	311	0	0	0	223	0	0	0	115	0	42	0	291	139	0	0	133	0	126	0	0	257	0	2237	0	156	0	0	0	0	0	0	0	0	0	227	269	111	2098	230	0	0	123	0	1841	220	170	0	0	0	95	0	113	0	91	0	0	268	177	83	72	67	0	481	0	179	0	1772	135	0	250	74	0	0	0	93	0	211	0	156	0	0	0	659	554	0	97	0	315	0	107	0	46	0	0	151	0	0	113	1967	0	0	169	0	624	0	89	0	434	0	59	0	0	0	163	108	67	0	0	77	0	161	0	0	0	0	0	0	76	82	0	339	0	0	0	0	487	140	0	0	0	0	324	155	0	0	0	0	0	96	0	0	0	0	0	0	0	0	0	0	0	276	0	0	0	0	0	0	0	135	0	0	0	0	0	0	0	0	0	0	0	0	106	37	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	104	312	0	0	43	0	0	280	0	0	0	0	72	0	0	0	91	0	106	24	0	0	0	0	166	0	0	0	0	0	0	0	80	0	0	0	0	0	0	0	39	0	0	0	130	0	0	0	180	0	92	0	0	0	0	64	0	63	62	0	51	0	0	0	0	0	0	0	0	95	0	0	0	144	0	0	0	0	0	0	0	0	0	0	0	0	154	0	13	65	0	0	0	0	0	0	0	0	0	0	86	0	117	373	0	0	140	0	0	134	0	88	0	0	0	0	51	0	0	0	84	0	0	123	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	83	0	0	0	0	0	0	52	0	0	0	0	136	0	0	0	0	0	0	0	0	88	0	0	0	8	0	0	0	26	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	4	0	0	10	0
HCB04-V2	2261	2120	3782	1573	465	503	3665	0	1219	342	0	2588	0	2120	285	702	0	887	0	636	0	215	0	502	174	198	212	0	0	0	414	564	0	144	563	0	225	0	807	0	0	180	177	319	0	89	206	0	0	0	679	282	0	0	84	149	1018	0	0	0	0	118	0	154	0	141	46	0	0	108	0	0	0	0	0	0	0	0	0	0	0	184	0	0	0	0	242	535	144	0	107	0	0	0	0	0	0	0	0	229	0	116	0	157	0	0	439	0	0	0	0	0	0	0	0	0	174	0	0	0	0	0	0	0	0	0	0	94	82	0	0	0	0	0	0	0	0	0	0	0	19	173	0	0	0	325	68	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	306	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	74	0	0	0	0	0	0	0	172	121	0	0	0	0	0	0	0	0	0	0	0	6	0	0	0	0	100	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	73	0	0	0	0	0	0	0	0	0	0	0	83	0	0	0	0	26	0	0	0	0	89	0	0	255	0	0	0	0	0	0	58	0	0	0	0	0	0	242	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	112	0	0	0	0	0	0	0	0	0	0	0	0	0	146	0	0	62	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	54	106	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	113	0	0	51	0	0	0	89	0	0	0	0	0	0	0	0	0	0	0	47	0	0	0	0	0	0	0	15	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	6	0	0	13	0	6	0	0	0	0

#taxonomy table

head(tax_table(psd5)) %>%
  kable(format = "html", col.names = colnames(tax_table(psd5))) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "320px")

	Kingdom	Phylum	Class	Order	Family	Genus
ASV1	Bacteria	Firmicutes	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus
ASV2	Bacteria	Firmicutes	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus
ASV3	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Blautia
ASV4	Bacteria	Firmicutes	Bacilli	Lactobacillales	Streptococcaceae	Streptococcus
ASV5	Bacteria	Proteobacteria	Gammaproteobacteria	Enterobacterales	Enterobacteriaceae	Escherichia-Shigella
ASV6	Bacteria	Proteobacteria	Gammaproteobacteria	Enterobacterales	Enterobacteriaceae	Escherichia-Shigella

#metadata table

as(sample_data(psd5), "data.frame") -> metad
 
metad %>%
  kable(format = "html", col.names = colnames(metad)) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "400px")

	SampleID	Edad_a_la_inclusion	Sexo	Respuesta_100	aloTPH_previo_si_no	Citogenética_alto_riesgo	Afect_Extramedular_PET_screening_si_no	Dias_producción	CRS_si_no	Grado_máx_CRS	HOSPITAL	Visit
CUN02-V2	CUN02-V2	43	MUJER	sCR	NO	NO	SI	120	SI	10	CUN	V2
CUN03-V2	CUN03-V2	45	HOMBRE	sCR	NO	SI	NO	100	SI	20	CUN	V2
HCB01-V2	HCB01-V2	66	HOMBRE	PR	NO	NO	SI	110	SI	20	HCB	V2
HCB02-V2	HCB02-V2	66	HOMBRE	PR	NO	SI	NO	110	SI	20	HCB	V2
HCB03-V2	HCB03-V2	60	HOMBRE	sCR	NO	SI	SI	110	SI	10	HCB	V2
HCB04-V2	HCB04-V2	55	MUJER	sCR	NO	NO	NO	110	SI	10	HCB	V2
HCB06-V2	HCB06-V2	73	MUJER	sCR	NO	NO	NO	90	SI	10	HCB	V2
HCB08-V2	HCB08-V2	62	MUJER	PR	NO	NO	SI	100	SI	10	HCB	V2
HCB09-V2	HCB09-V2	53	HOMBRE	sCR	NO	NO	NO	90	SI	10	HCB	V2
HCB10-V2	HCB10-V2	52	HOMBRE	sCR	SI	NO	SI	100	SI	10	HCB	V2
HCB12-V2	HCB12-V2	64	HOMBRE	sCR	NO	SI	NO	110	NO	NA	HCB	V2
HVA01-V2	HVA01-V2	58	HOMBRE	PR	NO	NO	NO	120	SI	20	HVA	V2
HVA02-V2	HVA02-V2	67	HOMBRE	sCR	NO	NO	SI	100	SI	20	HVA	V2
HVA03-V2	HVA03-V2	61	HOMBRE	sCR	SI	NO	SI	140	SI	10	HVA	V2
HVA04-V2	HVA04-V2	74	MUJER	sCR	NO	NO	NO	100	SI	10	HVA	V2
HVR01-V2	HVR01-V2	65	HOMBRE	sCR		NO	SI	110	SI	10	HVR	V2
HVR02-V2	HVR02-V2	46	HOMBRE	sCR		NO	SI	120	SI	10	HVR	V2
ISBS01-V2	ISBS01-V2	56	HOMBRE	sCR	NO	SI	SI	120	SI	10	ISBS	V2
ISBS03-V2	ISBS03-V2	70	MUJER	sCR	NO		SI	100	SI	10	ISBS	V2
ISBS06-V2	ISBS06-V2	60	HOMBRE	PR	NO	NO	NO	130	SI	10	ISBS	V2
ISBS08-V2	ISBS08-V2	53	HOMBRE	sCR	NO		SI	110	SI	20	ISBS	V2
ISBS09-V2	ISBS09-V2	65	MUJER	PR	NO	NO	SI	100	SI	10	ISBS	V2
ISBS10-V2	ISBS10-V2	68	HOMBRE	sCR	NO	SI	SI	140	NO	NA	ISBS	V2

#phylogenetic tree

plot_tree(psd5, method = "treeonly", ladderize = "left")

##Global summary

library(phyloseq)
library(microbiome)
summarize_phyloseq(psd5)

## Compositional = NO2

## 1] Min. number of reads = 146862] Max. number of reads = 1009443] Total number of reads = 7463814] Average number of reads = 32451.3478260875] Median number of reads = 329137] Sparsity = 0.7430333307314036] Any OTU sum to 1 or less? NO8] Number of singletons = 09] Percent of OTUs that are singletons 
##         (i.e. exactly one read detected across all samples)010] Number of sample variables are: 12SampleIDEdad_a_la_inclusionSexoRespuesta_100aloTPH_previo_si_noCitogenética_alto_riesgoAfect_Extramedular_PET_screening_si_noDias_producciónCRS_si_noGrado_máx_CRSHOSPITALVisit2

## [[1]]
## [1] "1] Min. number of reads = 14686"
## 
## [[2]]
## [1] "2] Max. number of reads = 100944"
## 
## [[3]]
## [1] "3] Total number of reads = 746381"
## 
## [[4]]
## [1] "4] Average number of reads = 32451.347826087"
## 
## [[5]]
## [1] "5] Median number of reads = 32913"
## 
## [[6]]
## [1] "7] Sparsity = 0.743033330731403"
## 
## [[7]]
## [1] "6] Any OTU sum to 1 or less? NO"
## 
## [[8]]
## [1] "8] Number of singletons = 0"
## 
## [[9]]
## [1] "9] Percent of OTUs that are singletons \n        (i.e. exactly one read detected across all samples)0"
## 
## [[10]]
## [1] "10] Number of sample variables are: 12"
## 
## [[11]]
##  [1] "SampleID"                              
##  [2] "Edad_a_la_inclusion"                   
##  [3] "Sexo"                                  
##  [4] "Respuesta_100"                         
##  [5] "aloTPH_previo_si_no"                   
##  [6] "Citogenética_alto_riesgo"              
##  [7] "Afect_Extramedular_PET_screening_si_no"
##  [8] "Dias_producción"                       
##  [9] "CRS_si_no"                             
## [10] "Grado_máx_CRS"                         
## [11] "HOSPITAL"                              
## [12] "Visit"

#group all tables

df <- psmelt(psd5)

head(df) %>%
  kable(format = "html", col.names = colnames(df)) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "400px")

	OTU	Sample	Abundance	SampleID	Edad_a_la_inclusion	Sexo	Respuesta_100	aloTPH_previo_si_no	Citogenética_alto_riesgo	Afect_Extramedular_PET_screening_si_no	Dias_producción	CRS_si_no	Grado_máx_CRS	HOSPITAL	Visit	Kingdom	Phylum	Class	Order	Family	Genus
4871	ASV29	ISBS03-V2	4933	ISBS03-V2	70	MUJER	sCR	NO		SI	100	SI	10	ISBS	V2	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Verrucomicrobiales	Akkermansiaceae	Akkermansia
7426	ASV39	ISBS03-V2	4171	ISBS03-V2	70	MUJER	sCR	NO		SI	100	SI	10	ISBS	V2	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Verrucomicrobiales	Akkermansiaceae	Akkermansia
7702	ASV40	ISBS03-V2	4020	ISBS03-V2	70	MUJER	sCR	NO		SI	100	SI	10	ISBS	V2	Bacteria	Verrucomicrobiota	Verrucomicrobiae	Verrucomicrobiales	Akkermansiaceae	Akkermansia
9466	ASV47	ISBS03-V2	3935	ISBS03-V2	70	MUJER	sCR	NO		SI	100	SI	10	ISBS	V2	Bacteria	Actinobacteriota	Actinobacteria	Bifidobacteriales	Bifidobacteriaceae	Bifidobacterium
5117	ASV3	HCB04-V2	3782	HCB04-V2	55	MUJER	sCR	NO	NO	NO	110	SI	10	HCB	V2	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Blautia
12070	ASV7	HCB04-V2	3665	HCB04-V2	55	MUJER	sCR	NO	NO	NO	110	SI	10	HCB	V2	Bacteria	Firmicutes	Clostridia	Lachnospirales	Lachnospiraceae	Blautia

##Distribution of the samples according to the metadata <> and <>

res <- plot_frequencies(sample_data(psd5), "Respuesta_100", "Sexo")
print(res$plot)

## <environment: 0x55d58a33a180>

# En formato tabla
res$data %>%
  kable(format = "html", col.names = colnames(res$data), digits = 2) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "300px")

Groups	Factor	n	pct	.group
PR	HOMBRE	4	66.67	1
PR	MUJER	2	33.33	1
sCR	HOMBRE	12	70.59	2
sCR	MUJER	5	29.41	2

#Taxa

taxa_sums(psd5)

##   ASV1   ASV2   ASV3   ASV4   ASV5   ASV6   ASV7   ASV8   ASV9  ASV10  ASV11 
##  16001  14310  12410  10325   8156   7768  10920   6885   7832   6216   6528 
##  ASV12  ASV13  ASV14  ASV15  ASV16  ASV17  ASV18  ASV19  ASV20  ASV21  ASV22 
##   8577   8130   7045   4841   7658   4986   5973   7514   6463   3963   4877 
##  ASV23  ASV24  ASV25  ASV26  ASV27  ASV28  ASV29  ASV30  ASV31  ASV32  ASV33 
##   5876   5391   4502   5060   3177   4933   6493   3813   4751   4545   3698 
##  ASV34  ASV35  ASV36  ASV37  ASV38  ASV39  ASV40  ASV41  ASV42  ASV43  ASV44 
##   4569   3529   3521   3332   5288   5702   5637   1072   3452   2834   3858 
##  ASV45  ASV46  ASV47  ASV48  ASV49  ASV50  ASV51  ASV52  ASV53  ASV54  ASV55 
##   4591   5149   4864   3673   2666   3375   2728   2836   3312    790   4453 
##  ASV56  ASV57  ASV58  ASV59  ASV60  ASV61  ASV62  ASV63  ASV64  ASV65  ASV66 
##   2917   2955   3666   3608   2976   3338   3337    895   3947   1999   2945 
##  ASV67  ASV68  ASV69  ASV70  ASV71  ASV72  ASV73  ASV74  ASV75  ASV76  ASV77 
##   3708    978   2811   2929   2802   2219   2466   2719   2789   2544   2847 
##  ASV78  ASV79  ASV80  ASV81  ASV82  ASV83  ASV84  ASV85  ASV86  ASV87  ASV88 
##   2556   1805   2649   2276   2034   2762   2174   3522   2308   2454   1718 
##  ASV89  ASV90  ASV91  ASV92  ASV93  ASV94  ASV95  ASV96  ASV97  ASV98  ASV99 
##   2236   1975   2689   2876   2240   2675   2609   2291   2035   2622   3095 
## ASV100 ASV101 ASV102 ASV103 ASV104 ASV105 ASV106 ASV107 ASV108 ASV109 ASV110 
##   2188   2879   2822   2208   1884   2089   2419   2361   2883   2181   2027 
## ASV111 ASV112 ASV113 ASV114 ASV115 ASV116 ASV117 ASV118 ASV119 ASV120 ASV121 
##   2262   2030   1472   2159   1794    604   1872   2710   1977   1938   2244 
## ASV122 ASV123 ASV124 ASV125 ASV126 ASV127 ASV128 ASV129 ASV130 ASV131 ASV132 
##   1437   1796   1399   1424   1923   1731   1890   1624   2010   2427   1658 
## ASV133 ASV134 ASV135 ASV136 ASV137 ASV138 ASV139 ASV140 ASV141 ASV142 ASV143 
##   1685   1927   1616   1365   2071   1738   1619   2087   1682   1331   1698 
## ASV144 ASV145 ASV146 ASV147 ASV148 ASV149 ASV150 ASV151 ASV152 ASV153 ASV154 
##   1830   1299   2257   1729   1157   1938   1722   1478   2136   1593   1551 
## ASV155 ASV156 ASV157 ASV158 ASV159 ASV160 ASV161 ASV162 ASV163 ASV164 ASV165 
##   1506   1932   1472   1405   1950   1312   1842   2161    245   1310   1885 
## ASV166 ASV167 ASV168 ASV169 ASV170 ASV171 ASV172 ASV173 ASV174 ASV175 ASV176 
##   1255   1699   1639    781    710   1575   1572   1849   1417   1260   1379 
## ASV177 ASV178 ASV179 ASV180 ASV181 ASV182 ASV183 ASV184 ASV185 ASV186 ASV187 
##   1244   1303   1896   1577   1249   1758    958   1137    951   1170   1079 
## ASV188 ASV189 ASV190 ASV191 ASV192 ASV193 ASV194 ASV195 ASV196 ASV197 ASV198 
##   1294   1326   1470   1223   1042   1733   1057   1014   1705   1231   1483 
## ASV199 ASV200 ASV201 ASV202 ASV203 ASV204 ASV205 ASV206 ASV207 ASV208 ASV209 
##    167   1158   1308   1371   1176   1113   1299    629   1045    628    678 
## ASV210 ASV211 ASV212 ASV213 ASV214 ASV215 ASV216 ASV217 ASV218 ASV219 ASV220 
##    615    697    981   1096    496   1072    923   1114    525    843   1174 
## ASV221 ASV222 ASV223 ASV224 ASV225 ASV226 ASV227 ASV228 ASV229 ASV230 ASV231 
##   1057    569   1058    809    435   1348    775   1176   1314    496    776 
## ASV232 ASV233 ASV234 ASV235 ASV236 ASV237 ASV238 ASV239 ASV240 ASV241 ASV242 
##   1282    901    494    700    769    874    407    712    912    826    562 
## ASV243 ASV244 ASV245 ASV246 ASV247 ASV248 ASV249 ASV250 ASV251 ASV252 ASV253 
##    484   1036    698    506    465   1179    778   1111    438    843   1004 
## ASV254 ASV255 ASV256 ASV257 ASV258 ASV259 ASV260 ASV261 ASV262 ASV263 ASV264 
##    645    760    636   1035    742    861    699    816    856    769    749 
## ASV265 ASV266 ASV267 ASV268 ASV269 ASV270 ASV271 ASV272 ASV273 ASV274 ASV275 
##    843    383    820    687    475    458    548   1043    569    611    499 
## ASV276 ASV277 ASV278 ASV279 ASV280 ASV281 ASV282 ASV283 ASV284 ASV285 ASV286 
##    436    648    988    919    728    365    604    448    547    654    938 
## ASV287 ASV288 ASV289 ASV290 ASV291 ASV292 ASV293 ASV294 ASV295 ASV296 ASV297 
##    542    579    651    526    418    360    538    382    913    326    709 
## ASV298 ASV299 ASV300 ASV301 ASV302 ASV303 ASV304 ASV305 ASV306 ASV307 ASV308 
##    264    355    656    686    544    582    753    501    689    579    406 
## ASV309 ASV310 ASV311 ASV312 ASV313 ASV314 ASV315 ASV316 ASV317 ASV318 ASV319 
##    492    354    725    372    661    317    574    303    470    683    682 
## ASV320 ASV321 ASV322 ASV323 ASV324 ASV325 ASV326 ASV327 ASV328 ASV329 ASV330 
##    650    423    540    719    453    408    603    471    422    142    686 
## ASV331 ASV332 ASV333 ASV334 ASV335 ASV336 ASV337 ASV338 ASV339 ASV340 ASV341 
##    225    305    456    647    499    424    270    247    599    424    466 
## ASV342 ASV343 ASV344 ASV345 ASV346 ASV347 ASV348 ASV349 ASV350 ASV351 ASV352 
##    445    463    325    562    585    334    335    295    559    426    660 
## ASV353 ASV354 ASV355 ASV356 ASV357 ASV358 ASV359 ASV360 ASV361 ASV362 ASV363 
##    507    536    562    583    144    633    455    540    425    223    462 
## ASV364 ASV365 ASV366 ASV367 ASV368 ASV369 ASV370 ASV371 ASV372 ASV373 ASV374 
##    347    501    536    406    405    511    314    415    301    395    341 
## ASV375 ASV376 ASV377 ASV378 ASV379 ASV380 ASV381 ASV382 ASV383 ASV384 ASV385 
##    158    269    404    226    252    396    543    457    531    385    301 
## ASV386 ASV387 ASV388 ASV389 ASV390 ASV391 ASV392 ASV393 ASV394 ASV395 ASV396 
##    521    329    385    345    503    502    413    493    343    357    488 
## ASV397 ASV398 ASV399 ASV400 ASV401 ASV402 ASV403 ASV404 ASV405 ASV406 ASV407 
##    137    139    296    305    401    394    462    318    306    279    233 
## ASV408 ASV409 ASV410 ASV411 ASV412 ASV413 ASV414 ASV415 ASV416 ASV417 ASV418 
##    330    227    192    444    281    436    350    282    430    340    330 
## ASV419 ASV420 ASV421 ASV422 ASV423 ASV424 ASV425 ASV426 ASV427 ASV428 ASV429 
##    427    297    358    336    256    307    161    229    349    325    193 
## ASV430 ASV431 ASV432 ASV433 ASV434 ASV435 ASV436 ASV437 ASV438 ASV439 ASV440 
##    380    273    344    392    245    284    277    163    356    260    248 
## ASV441 ASV442 ASV443 ASV444 ASV445 ASV446 ASV447 ASV448 ASV449 ASV450 ASV451 
##    369    259    187    284    279    359    359    250    356    355    126 
## ASV452 ASV453 ASV454 ASV455 ASV456 ASV457 ASV458 ASV459 ASV460 ASV461 ASV462 
##    349    179    347    224    208    249    220    215    195    263    215 
## ASV463 ASV464 ASV465 ASV466 ASV467 ASV468 ASV469 ASV470 ASV471 ASV472 ASV473 
##    229    330    245    245    299    197    323    265    264    213    200 
## ASV474 ASV475 ASV476 ASV477 ASV478 ASV479 ASV480 ASV481 ASV482 ASV483 ASV484 
##     85    176    253    182    308    246    300     99    131    294    262 
## ASV485 ASV486 ASV487 ASV488 ASV489 ASV490 ASV491 ASV492 ASV493 ASV494 ASV495 
##    286    284    284    278    246    185    272    159    266     80    263 
## ASV496 ASV497 ASV498 ASV499 ASV500 ASV501 ASV502 ASV503 ASV504 ASV505 ASV506 
##    178    257     90    250    205    143    188    245     91    124    159 
## ASV507 ASV508 ASV509 ASV510 ASV511 ASV512 ASV513 ASV514 ASV515 ASV516 ASV517 
##    236    201    233    169     88    228     71    221    179     18    175 
## ASV518 ASV519 ASV520 ASV521 ASV522 ASV523 ASV524 ASV525 ASV526 ASV527 ASV528 
##    206    129    201    138    178    193    132    189    106    188    188 
## ASV529 ASV530 ASV531 ASV532 ASV533 ASV534 ASV535 ASV536 ASV537 ASV538 ASV539 
##    156    162    147    154    154    147    120     65    128    121     99 
## ASV540 ASV541 ASV542 ASV543 ASV544 ASV545 ASV546 ASV547 ASV548 ASV549 ASV550 
##    108    104    103     68    100     26     90     55     75     80     65 
## ASV551 ASV552 ASV553 ASV554 ASV555 ASV556 ASV557 
##     54     50     42     36     35     34     23

kable(tax_table(psd5)[,6], format = "markdown")

	Genus
ASV1	Streptococcus
ASV2	Streptococcus
ASV3	Blautia
ASV4	Streptococcus
ASV5	Escherichia-Shigella
ASV6	Escherichia-Shigella
ASV7	Blautia
ASV8	Ruminococcus
ASV9	Streptococcus
ASV10	Escherichia-Shigella
ASV11	Ruminococcus
ASV12	Blautia
ASV13	Collinsella
ASV14	Blautia
ASV15	Escherichia-Shigella
ASV16	Bifidobacterium
ASV17	Ruminococcus
ASV18	Streptococcus
ASV19	Collinsella
ASV20	Bifidobacterium
ASV21	Ruminococcus
ASV22	Streptococcus
ASV23	Collinsella
ASV24	Bifidobacterium
ASV25	Streptococcus
ASV26	Anaerostipes
ASV27	Escherichia-Shigella
ASV28	Ruminococcus
ASV29	Akkermansia
ASV30	Subdoligranulum
ASV31	Collinsella
ASV32	Bifidobacterium
ASV33	Streptococcus
ASV34	Ruminococcus
ASV35	Streptococcus
ASV36	Streptococcus
ASV37	Subdoligranulum
ASV38	Streptococcus
ASV39	Akkermansia
ASV40	Akkermansia
ASV41	Streptococcus
ASV42	Blautia
ASV43	Ruminococcus
ASV44	Anaerostipes
ASV45	Streptococcus
ASV46	Akkermansia
ASV47	Bifidobacterium
ASV48	Ligilactobacillus
ASV49	Streptococcus
ASV50	Anaerostipes
ASV51	Subdoligranulum
ASV52	Blautia
ASV53	Ruminococcus
ASV54	Streptococcus
ASV55	Bifidobacterium
ASV56	Bacteroides
ASV57	Blautia
ASV58	Streptococcus
ASV59	Streptococcus
ASV60	Blautia
ASV61	Ligilactobacillus
ASV62	Holdemanella
ASV63	HT002
ASV64	Blautia
ASV65	Escherichia-Shigella
ASV66	Bacteroides
ASV67	Streptococcus
ASV68	HT002
ASV69	Anaerostipes
ASV70	Enterococcus
ASV71	Enterococcus
ASV72	Bacteroides
ASV73	[Ruminococcus] gnavus group
ASV74	[Eubacterium] hallii group
ASV75	[Eubacterium] hallii group
ASV76	Dorea
ASV77	Holdemanella
ASV78	Bifidobacterium
ASV79	Ruminococcus
ASV80	Ruminococcus
ASV81	Blautia
ASV82	Bacteroides
ASV83	Holdemanella
ASV84	Subdoligranulum
ASV85	Bifidobacterium
ASV86	Bacteroides
ASV87	Blautia
ASV88	Ruminococcus
ASV89	Dorea
ASV90	Erysipelatoclostridium
ASV91	NA
ASV92	Streptococcus
ASV93	Subdoligranulum
ASV94	[Ruminococcus] gnavus group
ASV95	Ligilactobacillus
ASV96	[Eubacterium] hallii group
ASV97	Blautia
ASV98	Agathobacter
ASV99	Blautia
ASV100	Subdoligranulum
ASV101	Bacteroides
ASV102	Bacteroides
ASV103	NA
ASV104	Erysipelatoclostridium
ASV105	Bifidobacterium
ASV106	Streptococcus
ASV107	Holdemanella
ASV108	Bifidobacterium
ASV109	Enterococcus
ASV110	Bifidobacterium
ASV111	NA
ASV112	Faecalibacterium
ASV113	Bacteroides
ASV114	Agathobacter
ASV115	Subdoligranulum
ASV116	HT002
ASV117	[Eubacterium] hallii group
ASV118	Akkermansia
ASV119	Agathobacter
ASV120	[Eubacterium] hallii group
ASV121	Bacteroides
ASV122	Erysipelatoclostridium
ASV123	Blautia
ASV124	Bifidobacterium
ASV125	Bifidobacterium
ASV126	Ruminococcus
ASV127	Dorea
ASV128	Ligilactobacillus
ASV129	NA
ASV130	Faecalibacterium
ASV131	Blautia
ASV132	Bacteroides
ASV133	Fusicatenibacter
ASV134	Streptococcus
ASV135	Romboutsia
ASV136	Subdoligranulum
ASV137	Dorea
ASV138	Enterococcus
ASV139	Subdoligranulum
ASV140	Dialister
ASV141	Faecalibacterium
ASV142	Bacteroides
ASV143	Bifidobacterium
ASV144	Faecalibacterium
ASV145	Bifidobacterium
ASV146	Akkermansia
ASV147	Streptococcus
ASV148	Erysipelatoclostridium
ASV149	Bacteroides
ASV150	Blautia
ASV151	NA
ASV152	Bacteroides
ASV153	[Ruminococcus] torques group
ASV154	NA
ASV155	[Ruminococcus] gnavus group
ASV156	Blautia
ASV157	Fusicatenibacter
ASV158	[Eubacterium] hallii group
ASV159	Blautia
ASV160	Blautia
ASV161	Bacteroides
ASV162	Bacteroides
ASV163	Lacticaseibacillus
ASV164	Bacteroides
ASV165	Bacteroides
ASV166	[Eubacterium] hallii group
ASV167	Dorea
ASV168	Agathobacter
ASV169	Lactobacillus
ASV170	Lactobacillus
ASV171	Streptococcus
ASV172	Blautia
ASV173	Bacteroides
ASV174	Blautia
ASV175	Dorea
ASV176	NA
ASV177	Fusicatenibacter
ASV178	NA
ASV179	Klebsiella
ASV180	Bacteroides
ASV181	[Ruminococcus] torques group
ASV182	Akkermansia
ASV183	Faecalibacterium
ASV184	Faecalibacterium
ASV185	Bifidobacterium
ASV186	Blautia
ASV187	Bacteroides
ASV188	Ligilactobacillus
ASV189	Blautia
ASV190	Dorea
ASV191	Coprococcus
ASV192	[Eubacterium] hallii group
ASV193	Klebsiella
ASV194	Romboutsia
ASV195	Faecalibacterium
ASV196	Akkermansia
ASV197	Ruminococcus
ASV198	Bacteroides
ASV199	Lacticaseibacillus
ASV200	[Ruminococcus] gnavus group
ASV201	[Clostridium] innocuum group
ASV202	Bacteroides
ASV203	Faecalibacterium
ASV204	Ruminococcus
ASV205	Streptococcus
ASV206	Alistipes
ASV207	Lachnoclostridium
ASV208	Limosilactobacillus
ASV209	Enterococcus
ASV210	Enterococcus
ASV211	Fusicatenibacter
ASV212	Coprococcus
ASV213	Blautia
ASV214	Alistipes
ASV215	Streptococcus
ASV216	[Clostridium] innocuum group
ASV217	Dorea
ASV218	Limosilactobacillus
ASV219	[Ruminococcus] torques group
ASV220	Anaerostipes
ASV221	Streptococcus
ASV222	NA
ASV223	Blautia
ASV224	Alistipes
ASV225	Lactobacillus
ASV226	Bacteroides
ASV227	Alistipes
ASV228	Bacteroides
ASV229	NA
ASV230	NA
ASV231	Faecalibacterium
ASV232	NA
ASV233	Ligilactobacillus
ASV234	Enterococcus
ASV235	Eggerthella
ASV236	Bacteroides
ASV237	[Ruminococcus] torques group
ASV238	Lactobacillus
ASV239	Romboutsia
ASV240	Blautia
ASV241	Streptococcus
ASV242	Faecalibacterium
ASV243	Collinsella
ASV244	UCG-002
ASV245	[Ruminococcus] torques group
ASV246	Erysipelatoclostridium
ASV247	Limosilactobacillus
ASV248	Klebsiella
ASV249	Alistipes
ASV250	Streptococcus
ASV251	Collinsella
ASV252	Streptococcus
ASV253	Blautia
ASV254	Romboutsia
ASV255	Streptococcus
ASV256	Parabacteroides
ASV257	Methanobrevibacter
ASV258	UBA1819
ASV259	[Ruminococcus] torques group
ASV260	Coprococcus
ASV261	UCG-002
ASV262	Blautia
ASV263	NA
ASV264	Hespellia
ASV265	Anaerostipes
ASV266	Enterococcus
ASV267	Roseburia
ASV268	Coprococcus
ASV269	Ruminococcus
ASV270	NA
ASV271	Streptococcus
ASV272	NA
ASV273	Blautia
ASV274	Monoglobus
ASV275	Sellimonas
ASV276	Faecalibacterium
ASV277	[Eubacterium] hallii group
ASV278	Streptococcus
ASV279	Senegalimassilia
ASV280	Bacteroides
ASV281	Collinsella
ASV282	Bifidobacterium
ASV283	Eggerthella
ASV284	Lachnoclostridium
ASV285	UBA1819
ASV286	Bacteroides
ASV287	Lachnoclostridium
ASV288	Parabacteroides
ASV289	Butyricicoccus
ASV290	Veillonella
ASV291	Erysipelatoclostridium
ASV292	Collinsella
ASV293	Erysipelatoclostridium
ASV294	NA
ASV295	Bifidobacterium
ASV296	Alistipes
ASV297	Bacteroides
ASV298	Limosilactobacillus
ASV299	Faecalibacterium
ASV300	Bifidobacterium
ASV301	Bacteroides
ASV302	[Clostridium] innocuum group
ASV303	Haemophilus
ASV304	Streptococcus
ASV305	Sellimonas
ASV306	Roseburia
ASV307	Alistipes
ASV308	Eggerthella
ASV309	Monoglobus
ASV310	Lachnoclostridium
ASV311	Bacteroides
ASV312	Lachnoclostridium
ASV313	Dorea
ASV314	NA
ASV315	Streptococcus
ASV316	NA
ASV317	NA
ASV318	UCG-002
ASV319	Hespellia
ASV320	Bacteroides
ASV321	[Eubacterium] hallii group
ASV322	Clostridium sensu stricto 1
ASV323	Methanobrevibacter
ASV324	Alistipes
ASV325	Anaerostipes
ASV326	Bacteroides
ASV327	Monoglobus
ASV328	Parabacteroides
ASV329	Subdoligranulum
ASV330	Senegalimassilia
ASV331	Erysipelatoclostridium
ASV332	NA
ASV333	Alistipes
ASV334	Coprococcus
ASV335	[Eubacterium] hallii group
ASV336	Streptococcus
ASV337	Veillonella
ASV338	NA
ASV339	Streptococcus
ASV340	Coprococcus
ASV341	Bifidobacterium
ASV342	Lachnoclostridium
ASV343	[Clostridium] innocuum group
ASV344	Lactococcus
ASV345	Bacteroides
ASV346	Blautia
ASV347	Bacteroides
ASV348	Lachnoclostridium
ASV349	NA
ASV350	UCG-002
ASV351	Veillonella
ASV352	Haemophilus
ASV353	Roseburia
ASV354	Coprococcus
ASV355	Methanobrevibacter
ASV356	Streptococcus
ASV357	Lactococcus
ASV358	Bacteroides
ASV359	Christensenellaceae R-7 group
ASV360	Roseburia
ASV361	Anaerostipes
ASV362	NA
ASV363	UBA1819
ASV364	Faecalibacterium
ASV365	Bacteroides
ASV366	Senegalimassilia
ASV367	[Ruminococcus] torques group
ASV368	[Eubacterium] hallii group
ASV369	Blautia
ASV370	Veillonella
ASV371	Bacteroides
ASV372	[Eubacterium] hallii group
ASV373	Bacteroides
ASV374	Alistipes
ASV375	Lachnoclostridium
ASV376	Family XIII AD3011 group
ASV377	Bacteroides
ASV378	[Eubacterium] hallii group
ASV379	Eggerthella
ASV380	Subdoligranulum
ASV381	Blautia
ASV382	Bacteroides
ASV383	Coprococcus
ASV384	NA
ASV385	Bifidobacterium
ASV386	[Ruminococcus] torques group
ASV387	Erysipelatoclostridium
ASV388	UCG-005
ASV389	Flavonifractor
ASV390	Bacteroides
ASV391	Bacteroides
ASV392	Subdoligranulum
ASV393	[Ruminococcus] gauvreauii group
ASV394	Christensenellaceae R-7 group
ASV395	NA
ASV396	Streptococcus
ASV397	Lachnoclostridium
ASV398	Sellimonas
ASV399	Anaerostipes
ASV400	Blautia
ASV401	Incertae Sedis
ASV402	Coprococcus
ASV403	Lachnospiraceae NK4A136 group
ASV404	Sellimonas
ASV405	Monoglobus
ASV406	Phascolarctobacterium
ASV407	Incertae Sedis
ASV408	Slackia
ASV409	[Eubacterium] brachy group
ASV410	Bacteroides
ASV411	Christensenellaceae R-7 group
ASV412	Agathobacter
ASV413	Blautia
ASV414	Senegalimassilia
ASV415	NA
ASV416	Hespellia
ASV417	Alistipes
ASV418	Streptococcus
ASV419	Eubacterium
ASV420	UCG-002
ASV421	NA
ASV422	Clostridium sensu stricto 1
ASV423	Collinsella
ASV424	Streptococcus
ASV425	NA
ASV426	Family XIII AD3011 group
ASV427	Bacteroides
ASV428	Butyricicoccus
ASV429	Erysipelatoclostridium
ASV430	Coprococcus
ASV431	NA
ASV432	Veillonella
ASV433	NA
ASV434	Family XIII AD3011 group
ASV435	Bacteroides
ASV436	UCG-005
ASV437	Incertae Sedis
ASV438	Streptococcus
ASV439	Blautia
ASV440	Bacteroides
ASV441	Incertae Sedis
ASV442	[Ruminococcus] gauvreauii group
ASV443	Eggerthella
ASV444	Bacteroides
ASV445	Clostridium sensu stricto 1
ASV446	Coprococcus
ASV447	UBA1819
ASV448	Bacteroides
ASV449	Blautia
ASV450	Blautia
ASV451	[Eubacterium] siraeum group
ASV452	Bacteroides
ASV453	Alistipes
ASV454	Christensenellaceae R-7 group
ASV455	Gemella
ASV456	Intestinibacter
ASV457	Incertae Sedis
ASV458	NA
ASV459	Sellimonas
ASV460	Streptococcus
ASV461	Roseburia
ASV462	Subdoligranulum
ASV463	Bifidobacterium
ASV464	Lachnoclostridium
ASV465	Parabacteroides
ASV466	Blautia
ASV467	Incertae Sedis
ASV468	Roseburia
ASV469	Alistipes
ASV470	Anaerostipes
ASV471	Slackia
ASV472	Christensenellaceae R-7 group
ASV473	NK4A214 group
ASV474	Coprobacillus
ASV475	NA
ASV476	Oscillibacter
ASV477	Faecalibacterium
ASV478	[Eubacterium] eligens group
ASV479	Streptococcus
ASV480	Intestinimonas
ASV481	Latilactobacillus
ASV482	NA
ASV483	Subdoligranulum
ASV484	Blautia
ASV485	Parabacteroides
ASV486	Christensenellaceae R-7 group
ASV487	NA
ASV488	Olsenella
ASV489	Oscillibacter
ASV490	NA
ASV491	Lachnoclostridium
ASV492	Roseburia
ASV493	Parabacteroides
ASV494	Latilactobacillus
ASV495	Family XIII AD3011 group
ASV496	NA
ASV497	Blautia
ASV498	Flavonifractor
ASV499	Bacteroides
ASV500	DTU089
ASV501	Odoribacter
ASV502	CAG-56
ASV503	Christensenellaceae R-7 group
ASV504	Blautia
ASV505	Streptococcus
ASV506	Roseburia
ASV507	UCG-002
ASV508	Incertae Sedis
ASV509	Terrisporobacter
ASV510	Oscillibacter
ASV511	Lachnoclostridium
ASV512	Lachnospiraceae NK4A136 group
ASV513	Lachnoclostridium
ASV514	Anaerostipes
ASV515	Incertae Sedis
ASV516	Solobacterium
ASV517	Alistipes
ASV518	[Ruminococcus] torques group
ASV519	Intestinimonas
ASV520	Christensenellaceae R-7 group
ASV521	[Eubacterium] siraeum group
ASV522	Family XIII AD3011 group
ASV523	Christensenellaceae R-7 group
ASV524	Streptococcus
ASV525	Bacteroides
ASV526	Family XIII AD3011 group
ASV527	Coprobacillus
ASV528	UBA1819
ASV529	Oscillibacter
ASV530	Christensenellaceae R-7 group
ASV531	NK4A214 group
ASV532	Incertae Sedis
ASV533	[Eubacterium] nodatum group
ASV534	UC5-1-2E3
ASV535	Christensenellaceae R-7 group
ASV536	Family XIII UCG-001
ASV537	NA
ASV538	Alistipes
ASV539	[Eubacterium] brachy group
ASV540	Lachnospira
ASV541	Holdemania
ASV542	Bilophila
ASV543	Scardovia
ASV544	Family XIII AD3011 group
ASV545	TM7x
ASV546	NA
ASV547	Anaerotruncus
ASV548	Actinomyces
ASV549	Solobacterium
ASV550	NA
ASV551	Lachnospiraceae FCS020 group
ASV552	Turicibacter
ASV553	Bifidobacterium
ASV554	NA
ASV555	UCG-005
ASV556	Parabacteroides
ASV557	Candidatus Soleaferrea

#ASVs of the same gender at the sequence level are different even though they correspond to the same gender. Visualization of sequences of the taxonomic range in phylogenetic trees

# cluster by gender
psd5.genus = tax_glom(psd5, "Genus", NArm = FALSE)
# Luego por altura en el árbol filogenético
h1 = 0.4
psd5.tip = tip_glom(psd5, h = h1)
# plot a comparison to visualize the differences
multiPlotTitleTextSize = 15
p2tree = plot_tree(psd5, method = "treeonly",
                   ladderize = "left",
                   title = "Sin aglomeración") +
  theme(plot.title = element_text(size = multiPlotTitleTextSize))
p3tree = plot_tree(psd5.genus, method = "treeonly",
                   ladderize = "left", title = "A nivel de género") +
  theme(plot.title = element_text(size = multiPlotTitleTextSize))
p4tree = plot_tree(psd5.tip, method = "treeonly",
                   ladderize = "left", title = "Por altura") +
  theme(plot.title = element_text(size = multiPlotTitleTextSize))

# plot all trees 
grid.arrange(nrow = 1, p2tree, p3tree, p4tree)

####DIVERSITY ANALYSIS####

#Richness Summary Graphics
plot_richness(psd5)

## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.

#Richness Summary Graphics group by <<Respuesta_100>>
(p = plot_richness(psd5, x = "Respuesta_100"))

## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.

#Richness Summary Graphics group by <<Respuesta_100>>
plot_richness(psd5, color = "Respuesta_100", x = "Respuesta_100", measures = c("Observed", "Chao1", "Shannon")) + geom_boxplot(aes(fill = Respuesta_100), alpha=.7) + scale_color_manual(values = c("#a6cee3", "#b2df8a")) + scale_fill_manual(values = c("#a6cee3", "#b2df8a"))

## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
## 
## We recommended that you find the un-trimmed data and retry.

#Check if there is any significant effect of alpha diversity based on response type

# save a dataframe with the alpha diversity measures
alpha_pd <- estimate_pd(psd5)

## Calculating Faiths PD-index...

# combine the metadata with alpha.diversity
data <- cbind(sample_data(psd5), alpha_pd) 

# calculate an ANOVA
psd5.anova <- aov(PD ~ Respuesta_100, data) 
# install.packages("xtable")
library(xtable)
psd5.anova.table <- xtable(psd5.anova)
kable(psd5.anova.table, caption = "Tabla ANOVA", digits = 5, format = "markdown")

Tabla ANOVA
	Df	Sum Sq	Mean Sq	F value	Pr(>F)
Respuesta_100	1	1.38598	1.38598	0.0554	0.81619
Residuals	21	525.32500	25.01548	NA	NA

#Mann-Whitney-Wilcoxon Test
psd5.wilcox <- wilcox.test(PD ~ Respuesta_100, data) 
psd5.wilcox

## 
##  Wilcoxon rank sum exact test
## 
## data:  PD by Respuesta_100
## W = 57, p-value = 0.7077
## alternative hypothesis: true location shift is not equal to 0

#Tablle of divsersity scores
tab <- global(psd5, index = "all")

## Warning: The microbiome::global function is deprecated. 
##         Use the function microbiome::alpha instead.

## Observed richness

## Other forms of richness

## Diversity

## Evenness

## Dominance

## Rarity

head(tab) %>%
  kable(format = "html", col.names = colnames(tab), digits = 2) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "300px")

	observed	chao1	diversity_inverse_simpson	diversity_gini_simpson	diversity_shannon	diversity_fisher	diversity_coverage	evenness_camargo	evenness_pielou	evenness_simpson	evenness_evar	evenness_bulla	dominance_dbp	dominance_dmn	dominance_absolute	dominance_relative	dominance_simpson	dominance_core_abundance	dominance_gini	rarity_log_modulo_skewness	rarity_low_abundance	rarity_rare_abundance
CUN02-V2	121	121	44.17	0.98	4.23	15.62	17	0.50	0.88	0.37	0.45	0.57	0.06	0.11	2081	0.06	0.02	0.07	0.90	1.59	0.03	0.24
CUN03-V2	188	188	76.19	0.99	4.82	27.05	37	0.50	0.92	0.41	0.59	0.68	0.05	0.10	1391	0.05	0.01	0.24	0.82	1.74	0.04	0.25
HCB01-V2	193	193	89.62	0.99	4.87	26.84	37	0.51	0.93	0.46	0.57	0.66	0.03	0.07	1192	0.03	0.01	0.16	0.81	1.65	0.06	0.31
HCB02-V2	88	88	29.97	0.97	3.94	12.28	13	0.77	0.88	0.34	0.55	0.61	0.10	0.19	1639	0.10	0.03	0.04	0.92	1.76	0.01	0.58
HCB03-V2	126	126	36.74	0.97	4.20	16.75	15	0.49	0.87	0.29	0.49	0.58	0.07	0.14	2237	0.07	0.03	0.21	0.90	1.67	0.02	0.19
HCB04-V2	83	83	22.14	0.95	3.61	10.14	8	0.33	0.82	0.27	0.32	0.49	0.10	0.20	3782	0.10	0.05	0.66	0.95	1.59	0.01	0.10

library(microbiome)
# Richness
tab <- richness(psd5)
# Dominance
tab <- dominance(psd5, index = "all")
# Rarity
tab <- rarity(psd5, index = "all")
# Coverage
tab <- coverage(psd5, threshold = 0.5)
# inequality
tab <- inequality(psd5)
# evennes
tab <- evenness(psd5, "all")

library(ggpubr)
# Generate a `phyloseq` object without taxa that adds 0 reads

psd5.2 <- prune_taxa(taxa_sums(psd5) > 0, psd5)
# calculate the diversity indices
tab <- diversities(psd5.2, index = "all")

## Warning: 'diversities' is deprecated.
## Use 'diversity' instead.
## See help("Deprecated") and help("The microbiome::diversities function has been 
##         replaced by function microbiome::alpha. Update your code accordingly.-deprecated").

## Observed richness

## Other forms of richness

## Diversity

## Evenness

## Dominance

## Rarity

# visualize the results table
head(tab) %>%
  kable(format = "html", col.names = colnames(tab), digits = 2) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "310px")

	observed	chao1	diversity_inverse_simpson	diversity_gini_simpson	diversity_shannon	diversity_fisher	diversity_coverage	evenness_camargo	evenness_pielou	evenness_simpson	evenness_evar	evenness_bulla	dominance_dbp	dominance_dmn	dominance_absolute	dominance_relative	dominance_simpson	dominance_core_abundance	dominance_gini	rarity_log_modulo_skewness	rarity_low_abundance	rarity_rare_abundance
CUN02-V2	121	121	44.17	0.98	4.23	15.62	17	0.50	0.88	0.37	0.45	0.57	0.06	0.11	2081	0.06	0.02	0.07	0.90	1.59	0.03	0.24
CUN03-V2	188	188	76.19	0.99	4.82	27.05	37	0.50	0.92	0.41	0.59	0.68	0.05	0.10	1391	0.05	0.01	0.24	0.82	1.74	0.04	0.25
HCB01-V2	193	193	89.62	0.99	4.87	26.84	37	0.51	0.93	0.46	0.57	0.66	0.03	0.07	1192	0.03	0.01	0.16	0.81	1.65	0.06	0.31
HCB02-V2	88	88	29.97	0.97	3.94	12.28	13	0.77	0.88	0.34	0.55	0.61	0.10	0.19	1639	0.10	0.03	0.04	0.92	1.76	0.01	0.58
HCB03-V2	126	126	36.74	0.97	4.20	16.75	15	0.49	0.87	0.29	0.49	0.58	0.07	0.14	2237	0.07	0.03	0.21	0.90	1.67	0.02	0.19
HCB04-V2	83	83	22.14	0.95	3.61	10.14	8	0.33	0.82	0.27	0.32	0.49	0.10	0.20	3782	0.10	0.05	0.66	0.95	1.59	0.01	0.10

#extract metadata
psd5.2.meta <- meta(psd5.2)

head(psd5.2.meta) %>%
  kable(format = "html", col.names = colnames(psd5.2.meta), digits = 2) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "500px")

	SampleID	Edad_a_la_inclusion	Sexo	Respuesta_100	aloTPH_previo_si_no	Citogenética_alto_riesgo	Afect_Extramedular_PET_screening_si_no	Dias_producción	CRS_si_no	Grado_máx_CRS	HOSPITAL	Visit
CUN02-V2	CUN02-V2	43	MUJER	sCR	NO	NO	SI	120	SI	10	CUN	V2
CUN03-V2	CUN03-V2	45	HOMBRE	sCR	NO	SI	NO	100	SI	20	CUN	V2
HCB01-V2	HCB01-V2	66	HOMBRE	PR	NO	NO	SI	110	SI	20	HCB	V2
HCB02-V2	HCB02-V2	66	HOMBRE	PR	NO	SI	NO	110	SI	20	HCB	V2
HCB03-V2	HCB03-V2	60	HOMBRE	sCR	NO	SI	SI	110	SI	10	HCB	V2
HCB04-V2	HCB04-V2	55	MUJER	sCR	NO	NO	NO	110	SI	10	HCB	V2

#add the diversity table to the metadata
psd5.2.meta$Shannon <- tab$diversity_shannon

# get the variables from our `phyloseq` object
rpp <- levels(psd5.2.meta$Respuesta_100)
#Create a list of what we want to compare

pares.rpp <- combn(seq_along(rpp), 2, simplify = FALSE, FUN = function(i)rpp[i])
# Print results
print(pares.rpp)

## [[1]]
## [1] "PR"  "sCR"

#Plot 
p1 <- ggviolin(psd5.2.meta, x = "Respuesta_100", y = "Shannon",
 add = "boxplot", fill = "Respuesta_100", palette = c("#a6cee3", "#b2df8a"))  
print(p1)

#Plot with compare means
p1 <- p1 + stat_compare_means(comparisons = pares.rpp)
print(p1)

###Beta diversity

#UniFrac
psd5.mds.unifrac <- ordinate(psd5, method = "MDS", distance = "unifrac")

## Warning in UniFrac(physeq, ...): Randomly assigning root as -- ASV415 -- in the
## phylogenetic tree in the data you provided.

evals <- psd5.mds.unifrac$values$Eigenvalues
pord1 <- plot_ordination(psd5, psd5.mds.unifrac, color = "Respuesta_100") +
  labs(col = "Respuesta_100") +
  coord_fixed(sqrt(evals[2] / evals[1]))

#Bray-Curtis
psd5.mds.bray <- ordinate(psd5, method = "MDS", distance = "bray")
evals <- psd5.mds.bray$values$Eigenvalues
pord2 <- plot_ordination(psd5, psd5.mds.bray, color = "Respuesta_100") +
  labs(col = "Respuesta_100") +
  coord_fixed(sqrt(evals[2] / evals[1]))

grid.arrange(pord1, pord2)

#DPCoA
psd5.dpcoa.unifrac <- ordinate(psd5, method = "DPCoA", distance = "dpcoa")
evals <- psd5.dpcoa.unifrac$eig
pord3 <- plot_ordination(psd5, psd5.dpcoa.unifrac, color = "Respuesta_100", shape = "Sexo") +
  labs(col = "Respuesta_100") +
  coord_fixed(sqrt(evals[2] / evals[1])) +
  scale_color_manual(values = c("#a6cee3", "#b2df8a")) + 
  scale_fill_manual(values = c("#a6cee3",  "#b2df8a")) +
  geom_point(size=4)

## Warning: non-unique values when setting 'row.names':

## Species coordinates not found directly in ordination object. Attempting weighted average (`vegan::wascores`)

pord3

#multidimensional scaling

library(tsnemicrobiota)
tsne_res <- tsne_phyloseq(psd5, distance= "dpcoa", perplexity = 8, verbose=0, rng_seed = 3901)

# Plot
pord4 <- plot_tsne_phyloseq(psd5, tsne_res, color = "Respuesta_100", shape = "Sexo") +
  geom_point(size=4) +
  scale_color_manual(values = c("#a6cee3",  "#b2df8a")) + 
  scale_fill_manual(values = c("#a6cee3",  "#b2df8a"))

grid.arrange(pord3, pord4)

## Warning: Use of `plot_df[[axes[1]]]` is discouraged. Use `.data[[axes[1]]]`
## instead.

## Warning: Use of `plot_df[[axes[2]]]` is discouraged. Use `.data[[axes[2]]]`
## instead.

## Warning: Use of `plot_df[[axes[1]]]` is discouraged. Use `.data[[axes[1]]]`
## instead.

## Warning: Use of `plot_df[[axes[2]]]` is discouraged. Use `.data[[axes[2]]]`
## instead.

###Abundance analysis

# Necesitamos obtener las taxa más abundantes, en este caso el top 15
top15 <- get_top_taxa(physeq_obj = psd5, n = 15, relative = T,
                       discard_other = T, other_label = "Other")
# Ya que no todas las taxa fueron clasificadas a nivel de especie, generamos etiquetas compuestas de distintos rangos taxonómicos para el gráfico
top15 <- name_taxa(top15, label = "", species = F, other_label = "Other")
# Finalmente graficamos
fantaxtic_bar(top15, color_by = "Family", label_by = "Genus", facet_by = NULL, grid_by = NULL, other_color = "Grey") -> ptop15

##                Level N.color.shades Central.color
## 1   Streptococcaceae              4       #6495ed
## 2    Lachnospiraceae              4       #ff7256
## 3 Enterobacteriaceae              2       #edbc64
## 4  Coriobacteriaceae              3       #8470ff
## 5 Bifidobacteriaceae              2       #8ee5ee

#Barplot of abundance
ptop15

#Groupo by Respuesta_100
fantaxtic_bar(top15, color_by = "Family", label_by = "Genus", facet_by = "Respuesta_100", grid_by = NULL, other_color = "Grey") -> ptop15.2

##                Level N.color.shades Central.color
## 1   Streptococcaceae              4       #6495ed
## 2    Lachnospiraceae              4       #ff7256
## 3 Enterobacteriaceae              2       #edbc64
## 4  Coriobacteriaceae              3       #8470ff
## 5 Bifidobacteriaceae              2       #8ee5ee

#Plot 
ptop15.2

library(ampvis2)
# extract the read counts table and the taxonomy table from the psd5 object
#  generate a copy so as not to overwrite psd5
obj <- psd5
# change the orientation of the otu_table
t(otu_table(obj)) -> otu_table(obj)
# Extraemos las tablas
otutable <- data.frame(OTU = rownames(phyloseq::otu_table(obj)@.Data),
                       phyloseq::otu_table(obj)@.Data,
                       phyloseq::tax_table(obj)@.Data,
                       check.names = FALSE
)
# extract la metadada
metadata <- data.frame(phyloseq::sample_data(obj), 
                       check.names = FALSE
)

#ampvis2 requires that 1) the taxonomic ranges are seven and go from Kingdom to Species and 2) the first column of the metadata is the identifier of each sample
# Then we duplicate the Genus column and rename it to Species
otutable$Species = otutable$Genus
# Reordenamos la metadata
metadata <- metadata[,c("SampleID", "Edad_a_la_inclusion", "Sexo",  "Respuesta_100",    "aloTPH_previo_si_no",  "Citogenética_alto_riesgo", "Afect_Extramedular_PET_screening_si_no",   "Dias_producción",  "CRS_si_no", "Grado_máx_CRS",   "HOSPITAL", "Visit")]

# generate the ampvis object

av2 <- amp_load(otutable, metadata)

#estructura de las comunidades en las curvas de abundancia
amp_rankabundance(av2, log10_x = T, group_by = "Sexo")

#estructura de las comunidades en las curvas de abundancia
amp_rankabundance(av2, log10_x = T, group_by = "Respuesta_100")

#estructura de las comunidades en las curvas de abundancia
amp_rankabundance(av2, log10_x = T, group_by = "HOSPITAL")

#heatmap of abundance
amp_heatmap(av2, 
            group_by = "Respuesta_100", 
            facet_by = "Sexo", 
            plot_values = TRUE,
            tax_show = 20,
            tax_aggregate = "Genus",
            tax_add = "Phylum",
            plot_colorscale = "sqrt",
            plot_legendbreaks = c(1, 5, 10)
            )

#Box Plots.
amp_boxplot(av2,
            group_by = "Respuesta_100",
            tax_show = 20,
            tax_aggregate = "Genus",
            tax_add = "Phylum",
            adjust_zero = T,
            plot_log = T) +
  scale_color_manual(values = c("#a6cee3", "#b2df8a")) + 
  scale_fill_manual(values = c("#a6cee3", "#b2df8a"))

#Venn Diagram
# some of these microorganisms are shared between all samples
amp_venn(av2, group_by = "Respuesta_100", cut_a = 0, cut_f = 50, text_size = 3)

#Differential abundance analysis

#we want to determine exactly which taxa is more represented in one condition versus another and to what extent. DESeq2 package,

# Create a DESeq2 object with the `phyloseq_to_deseq2` function
diagdds = phyloseq_to_deseq2(psd4, ~Respuesta_100)

## converting counts to integer mode

# We calculate the size factors as part of the normalization of the samples
# calculate geometric means prior to estimate size factors
gm_mean = function(x, na.rm=TRUE){
    exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
geoMeans = apply(counts(diagdds), 1, gm_mean)
diagdds = estimateSizeFactors(diagdds, geoMeans = geoMeans)

# We normalize and perform the parametric Wald test to determine differentially abundant taxa.
diagdds = DESeq(diagdds, test="Wald", fitType="local")

## using pre-existing size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 297 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

# Create a DESeq2 object with the `phyloseq_to_deseq2` function
diagdds = phyloseq_to_deseq2(psd5, ~Respuesta_100)

## converting counts to integer mode

# We calculate the size factors as part of the normalization of the samples
# calculate geometric means prior to estimate size factors
gm_mean = function(x, na.rm=TRUE){
    exp(sum(log(x[x > 0]), na.rm=na.rm) / length(x))
}
geoMeans = apply(counts(diagdds), 1, gm_mean)
diagdds = estimateSizeFactors(diagdds, geoMeans = geoMeans)

# We normalize and perform the parametric Wald test to determine differentially abundant taxa.
diagdds = DESeq(diagdds, test="Wald", fitType="local")

## using pre-existing size factors

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## fitting model and testing

## -- replacing outliers and refitting for 297 genes
## -- DESeq argument 'minReplicatesForReplace' = 7 
## -- original counts are preserved in counts(dds)

## estimating dispersions

## fitting model and testing

## estimating size factors
plotMA(diagdds)

## estimating dispersions
plotDispEsts(diagdds)

# save the results 
res = results(diagdds, cooksCutoff = FALSE)
# reorder 
res = res[order(res$padj, na.last=NA), ]

#especify contrast
res.PR.sCR <- results(diagdds, contrast=c("Respuesta_100","PR","sCR"))

#Visualizate generated table
res.PR.sCR  %>%
  kable(format = "html", col.names = colnames(res.PR.sCR )) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "400px")

	baseMean	log2FoldChange	lfcSE	stat	pvalue	padj
ASV1	729.6119002	-1.4986609	1.481991	-1.0112485	0.3118975	0.9759373
ASV2	647.9485881	-1.7206118	1.752171	-0.9819884	0.3261055	0.9919616
ASV3	567.9817458	-1.6673863	1.140755	-1.4616510	0.1438369	0.8239591
ASV4	469.9120778	-1.8321980	1.676678	-1.0927550	0.2745014	0.9181599
ASV5	335.6251807	0.1718067	1.840638	0.0933409	0.9256328	0.9919616
ASV6	318.8184377	0.2142279	1.845320	0.1160926	0.9075792	0.9919616
ASV7	502.3517870	-1.6004280	1.459003	-1.0969324	0.2726710	0.9181599
ASV8	231.8964389	0.1585241	3.285750	0.0482459	0.9615202	0.9919616
ASV9	354.9656963	-1.5880479	1.624910	-0.9773146	0.3284134	0.9919616
ASV10	256.4548414	0.2538495	1.762246	0.1440489	0.8854619	0.9919616
ASV11	220.4540728	0.2587368	3.285756	0.0787450	0.9372355	0.9919616
ASV12	394.4344611	-1.4672698	1.384096	-1.0600924	0.2891026	0.9506084
ASV13	387.6501461	0.1555353	1.799128	0.0864504	0.9311084	0.9919616
ASV14	324.1069347	-1.1306624	1.612926	-0.7010006	0.4833026	0.9919616
ASV15	200.5535457	0.0647585	1.799807	0.0359808	0.9712977	0.9919616
ASV16	344.9188676	-0.5298251	1.510569	-0.3507453	0.7257794	0.9919616
ASV17	98.1964810	1.9335166	3.286040	0.5884032	0.5562617	0.9919616
ASV18	276.9018480	-1.8157152	1.969764	-0.9217930	0.3566366	0.9919616
ASV19	356.5740209	0.3518504	1.923800	0.1828934	0.8548817	0.9919616
ASV20	287.9218237	-0.3332220	1.845926	-0.1805175	0.8567463	0.9919616
ASV21	79.4269991	2.0999301	3.286195	0.6390157	0.5228127	0.9919616
ASV22	269.4452622	2.7801277	3.056320	0.9096325	0.3630164	0.9919616
ASV23	278.4514632	0.3426924	1.952610	0.1755047	0.8606830	0.9919616
ASV24	243.9912163	-0.4321198	1.982185	-0.2180017	0.8274278	0.9919616
ASV25	203.9225727	4.3251729	3.203453	1.3501595	0.1769648	0.8239591
ASV26	212.9596269	-0.7757047	1.414269	-0.5484845	0.5833593	0.9919616
ASV27	130.9438077	0.2653469	2.032955	0.1305228	0.8961528	0.9919616
ASV28	220.1562885	1.1202239	3.120900	0.3589426	0.7196380	0.9919616
ASV29	77.4021602	2.0629353	3.286211	0.6277551	0.5301644	0.9919616
ASV30	160.9733415	0.9049913	2.271571	0.3983988	0.6903363	0.9919616
ASV31	224.1176959	0.2021397	1.976725	0.1022599	0.9185504	0.9919616
ASV32	205.1326821	-0.3727775	1.886760	-0.1975755	0.8433772	0.9919616
ASV33	171.2487384	-1.9942607	2.121064	-0.9402172	0.3471061	0.9919616
ASV34	203.2785965	0.8817885	3.285741	0.2683682	0.7884159	0.9919616
ASV35	170.9246972	3.9919722	2.560842	1.5588516	NA	NA
ASV36	161.5216435	-2.0737110	2.122820	-0.9768663	0.3286354	0.9919616
ASV37	142.5629474	0.7749435	2.295679	0.3375662	0.7356901	0.9919616
ASV38	224.9469085	-1.7253151	2.099813	-0.8216518	0.4112751	0.9919616
ASV39	76.5044086	1.7966353	3.175117	0.5658484	0.5714968	0.9919616
ASV40	80.0219857	2.1775840	3.286198	0.6626453	0.5075577	0.9919616
ASV41	13.1607033	1.2184245	3.289774	0.3703672	NA	NA
ASV42	164.2897896	-0.1230981	1.760244	-0.0699324	0.9442475	0.9919616
ASV43	96.3928250	0.3648341	3.286129	0.1110225	0.9115985	0.9919616
ASV44	154.1188333	-0.3650474	2.009368	-0.1816727	0.8558396	0.9919616
ASV45	191.8091766	-1.2744567	1.998382	-0.6377442	0.5236402	0.9919616
ASV46	74.9889042	2.1514907	3.286245	0.6546957	0.5126637	0.9919616
ASV47	52.8249013	0.9780768	3.286553	0.2975996	NA	NA
ASV48	21.0604389	-21.6189123	3.390043	-6.3771802	0.0000000	0.0000000
ASV49	146.8648450	3.0450581	2.889235	1.0539323	0.2919139	0.9517866
ASV50	136.5104373	-0.3760083	2.153932	-0.1745683	0.8614189	0.9919616
ASV51	117.6774520	0.6283688	2.589345	0.2426748	0.8082573	0.9919616
ASV52	133.4977730	-0.0250548	2.019118	-0.0124088	0.9900995	0.9952295
ASV53	147.0733955	1.3014031	3.285837	0.3960644	0.6920575	0.9919616
ASV54	4.6224055	-5.1142311	3.392380	-1.5075643	0.1316661	0.8239591
ASV55	52.8793226	1.1494281	3.286533	0.3497388	NA	NA
ASV56	104.9625328	-0.3409819	2.289277	-0.1489474	0.8815951	0.9919616
ASV57	108.5927540	-0.1162813	2.274199	-0.0511307	0.9592214	0.9919616
ASV58	73.4771208	-2.5146396	2.530601	-0.9936925	NA	NA
ASV59	154.7492050	-1.6283106	2.134479	-0.7628610	0.4455463	0.9919616
ASV60	136.6390867	-1.1831715	2.025412	-0.5841634	0.5591104	0.9919616
ASV61	20.9664993	-21.6141366	3.390046	-6.3757658	0.0000000	0.0000000
ASV62	55.9784879	-1.6604781	3.288210	-0.5049793	NA	NA
ASV63	26.8069736	25.5388291	3.323020	7.6854285	NA	NA
ASV64	161.7641332	1.5093733	2.403896	0.6278862	0.5300785	0.9919616
ASV65	81.6119418	0.1705516	2.576779	0.0661879	0.9472282	0.9919616
ASV66	105.2490022	0.1401011	2.433482	0.0575723	0.9540893	0.9919616
ASV67	74.9196349	-2.8076000	2.369384	-1.1849491	NA	NA
ASV68	26.8957111	5.7690474	3.025839	1.9065941	NA	NA
ASV69	117.3571525	-0.1675072	2.155159	-0.0777238	0.9380477	0.9919616
ASV70	12.5135518	2.4985861	3.290537	0.7593247	NA	NA
ASV71	13.4673147	4.1309178	3.294256	1.2539759	NA	NA
ASV72	92.3464960	1.0178174	2.001669	0.5084845	0.6111136	0.9919616
ASV73	101.0965922	0.8143772	2.287125	0.3560702	0.7217880	0.9919616
ASV74	129.7373560	-0.3492474	1.747405	-0.1998664	0.8415851	0.9919616
ASV75	134.3798926	-0.4867884	1.749176	-0.2782959	0.7807852	0.9919616
ASV76	122.7898792	-2.0713988	1.992039	-1.0398386	0.2984149	0.9568046
ASV77	37.8789280	-0.8495222	3.288146	-0.2583591	NA	NA
ASV78	138.2753466	-24.2427804	3.338224	-7.2621789	0.0000000	0.0000000
ASV79	34.5747689	2.0240758	3.287141	0.6157557	0.5380558	0.9919616
ASV80	117.7473657	1.3424219	3.285932	0.4085361	0.6828801	0.9919616
ASV81	105.9333316	0.2802472	1.885028	0.1486700	0.8818140	0.9919616
ASV82	82.2507713	1.2290321	2.395312	0.5130989	0.6078822	0.9919616
ASV83	46.6562218	-1.9613577	3.289373	-0.5962709	NA	NA
ASV84	93.4199713	0.6856970	2.585498	0.2652089	0.7908485	0.9919616
ASV85	36.6760082	1.6387017	3.286995	0.4985410	NA	NA
ASV86	86.0773333	0.2100168	2.419947	0.0867857	0.9308419	0.9919616
ASV87	90.7172497	-0.8728382	2.268956	-0.3846872	0.7004692	0.9919616
ASV88	32.9963492	1.9926359	3.287215	0.6061776	0.5443968	0.9919616
ASV89	102.9145267	-2.3340245	2.255945	-1.0346103	0.3008509	0.9568046
ASV90	78.3930056	0.7356288	1.958344	0.3756383	0.7071858	0.9919616
ASV91	111.0160229	0.9400974	3.285978	0.2860937	NA	NA
ASV92	123.6839801	-1.9560215	2.421463	-0.8077851	0.4192143	0.9919616
ASV93	108.1357923	-1.7595728	2.262041	-0.7778695	0.4366460	0.9919616
ASV94	102.5765207	0.5920580	2.015561	0.2937435	0.7689539	0.9919616
ASV95	18.9632444	-21.4660581	3.390115	-6.3319549	0.0000000	0.0000000
ASV96	109.5682686	-0.4250654	1.749458	-0.2429698	0.8080288	0.9919616
ASV97	94.5684637	0.1747212	1.759452	0.0993043	0.9208966	0.9919616
ASV98	130.1495422	-1.3724365	2.281659	-0.6015083	0.5475015	0.9919616
ASV99	117.7058418	1.6244231	3.285932	0.4943569	0.6210542	0.9919616
ASV100	102.9937411	-2.1619775	1.962247	-1.1017869	0.2705543	0.9181599
ASV101	37.4555346	0.1752158	2.825810	0.0620055	NA	NA
ASV102	30.3563940	0.8641839	3.236315	0.2670272	0.7894482	0.9919616
ASV103	104.7422520	-0.0628827	2.145861	-0.0293042	0.9766220	0.9919616
ASV104	77.2801867	1.0011760	1.953669	0.5124594	0.6083295	0.9919616
ASV105	112.9382241	-23.9538555	3.303113	-7.2519033	0.0000000	0.0000000
ASV106	116.4831763	-3.9291606	2.392126	-1.6425394	0.1004783	0.8011152
ASV107	29.7158086	-0.7888388	3.254288	-0.2423998	NA	NA
ASV108	27.2525697	1.3666612	3.287526	0.4157112	NA	NA
ASV109	8.5787506	2.9838815	3.293954	0.9058663	NA	NA
ASV110	108.7714319	-23.9037073	3.290335	-7.2648249	0.0000000	0.0000000
ASV111	50.5524922	3.6920541	3.287360	1.1231060	NA	NA
ASV112	89.3307534	-1.2525755	3.199650	-0.3914726	0.6954479	0.9919616
ASV113	61.5370469	1.2055573	2.538284	0.4749498	0.6348228	0.9919616
ASV114	107.2245667	-1.3185745	2.437291	-0.5410000	0.5885076	0.9919616
ASV115	83.7157143	-1.8099497	2.096484	-0.8633262	0.3879581	0.9919616
ASV116	15.0144655	24.7446523	3.323423	7.4455326	NA	NA
ASV117	89.4470181	-0.1933697	1.740902	-0.1110744	0.9115573	0.9919616
ASV118	36.4658140	2.0453866	3.287059	0.6222544	0.5337746	0.9919616
ASV119	96.7864149	-1.3176077	2.586825	-0.5093533	0.6105046	0.9919616
ASV120	91.3182811	0.5654276	1.983503	0.2850651	0.7755942	0.9919616
ASV121	29.7605845	0.5276219	2.813538	0.1875297	0.8512454	0.9919616
ASV122	59.7716010	0.7052857	2.086231	0.3380669	0.7353127	0.9919616
ASV123	70.0608725	-0.2748272	2.397299	-0.1146403	0.9087302	0.9919616
ASV124	37.3398965	1.9666688	3.287006	0.5983162	NA	NA
ASV125	72.4822695	0.0395099	3.286442	0.0120221	NA	NA
ASV126	85.8037812	0.8863308	3.286138	0.2697181	0.7873771	0.9919616
ASV127	81.2073174	-3.2216932	2.207768	-1.4592536	0.1444953	0.8239591
ASV128	14.0029345	-21.0527886	3.390374	-6.2095766	0.0000000	0.0000000
ASV129	75.4741929	0.5466655	2.408316	0.2269907	0.8204310	0.9919616
ASV130	84.6584330	-1.2258763	2.825675	-0.4338349	0.6644083	0.9919616
ASV131	96.1575862	1.2092372	3.222729	0.3752215	0.7074957	0.9919616
ASV132	60.4910383	-0.1026891	2.721032	-0.0377390	0.9698958	0.9919616
ASV133	77.1898867	-0.3008719	1.761643	-0.1707904	0.8643885	0.9919616
ASV134	91.4972903	-1.7801029	2.577559	-0.6906158	0.4898070	0.9919616
ASV135	63.6047328	0.7157825	1.963856	0.3644782	0.7155010	0.9919616
ASV136	57.8012037	0.8251714	2.565021	0.3217017	0.7476787	0.9919616
ASV137	102.0358370	-2.9711336	2.241957	-1.3252414	0.1850911	0.8254637
ASV138	6.5626060	2.0304048	3.294337	0.6163318	NA	NA
ASV139	76.8622643	-1.5698741	2.237371	-0.7016602	0.4828911	0.9919616
ASV140	6.2522866	-5.5479189	3.391603	-1.6357809	0.1018855	0.8011152
ASV141	81.6726360	0.5822274	2.425495	0.2400448	0.8102955	0.9919616
ASV142	55.4562266	1.1084292	2.386938	0.4643729	0.6423806	0.9919616
ASV143	90.5791912	-23.6475421	2.798983	-8.4486197	0.0000000	0.0000000
ASV144	82.3717268	-1.4105106	2.400063	-0.5876973	0.5567355	0.9919616
ASV145	35.1470250	1.7299924	3.287069	0.5263024	NA	NA
ASV146	33.1813171	2.1875452	3.201299	0.6833304	0.4943981	0.9919616
ASV147	76.9935825	-9.1714264	2.113207	-4.3400514	0.0000142	0.0003251
ASV148	47.4419723	0.4534362	2.245278	0.2019510	0.8399550	0.9919616
ASV149	26.7691252	-0.3444786	2.826824	-0.1218606	0.9030094	0.9919616
ASV150	68.3408444	0.5461074	2.555630	0.2136880	0.8307904	0.9919616
ASV151	69.3597547	0.6442583	2.733715	0.2356713	0.8136877	0.9919616
ASV152	22.9817179	0.7227008	3.288079	0.2197942	NA	NA
ASV153	67.3421817	-1.4197520	2.735684	-0.5189751	0.6037781	0.9919616
ASV154	33.1585356	3.3676116	3.287988	1.0242165	NA	NA
ASV155	58.9457682	0.8087817	2.405926	0.3361623	0.7367485	0.9919616
ASV156	80.0176808	2.9913610	2.494061	1.1993936	0.2303749	0.8845429
ASV157	68.6227101	0.0878385	2.404225	0.0365351	0.9708557	0.9919616
ASV158	62.2145850	0.6289519	2.238705	0.2809445	0.7787529	0.9919616
ASV159	77.5779727	1.0698959	3.173304	0.3371552	0.7359999	0.9919616
ASV160	48.8186867	-1.4597874	2.391080	-0.6105140	0.5415214	0.9919616
ASV161	23.9642660	0.6512551	3.288006	0.1980699	NA	NA
ASV162	17.1081541	1.0647190	3.288816	0.3237394	NA	NA
ASV163	12.4978267	2.3706809	3.290406	0.7204829	0.4712278	0.9919616
ASV164	46.7524583	0.2694123	3.286885	0.0819659	0.9346739	0.9919616
ASV165	26.2865534	-20.7220813	3.389912	-6.1128676	0.0000000	0.0000000
ASV166	57.2151267	0.5838973	2.096336	0.2785324	0.7806037	0.9919616
ASV167	82.8643127	-2.9766661	2.534423	-1.1744947	0.2401969	0.8845429
ASV168	79.9113593	-1.1468395	2.282712	-0.5024022	0.6153846	0.9919616
ASV169	30.5030962	2.4961217	2.952891	0.8453145	0.3979353	0.9919616
ASV170	28.6949520	2.4676398	3.287656	0.7505773	0.4529071	0.9919616
ASV171	65.7534356	-8.9438292	2.569298	-3.4810399	0.0004995	0.0107664
ASV172	65.3240072	3.8085358	3.287009	1.1586630	0.2465936	0.8941898
ASV173	27.8082679	-21.4318902	3.389883	-6.3223099	0.0000000	0.0000000
ASV174	57.5447541	0.6523279	2.555788	0.2552355	0.7985412	0.9919616
ASV175	57.5698035	-1.9116103	2.226085	-0.8587318	0.3904885	0.9919616
ASV176	31.2072349	2.5723920	3.207363	0.8020271	NA	NA
ASV177	58.1177857	0.0114375	2.261235	0.0050581	0.9959643	0.9968603
ASV178	61.4946545	0.2436578	2.571961	0.0947362	0.9245244	0.9919616
ASV179	50.0383825	26.3335235	3.322781	7.9251453	NA	NA
ASV180	19.9949146	-0.6167675	3.025000	-0.2038901	NA	NA
ASV181	21.3556423	1.6635071	3.288104	0.5059168	0.6129151	0.9919616
ASV182	23.8275322	1.9593435	3.241651	0.6044277	0.5455594	0.9919616
ASV183	40.2638844	-1.9626797	2.740299	-0.7162284	0.4738504	0.9919616
ASV184	53.7619259	1.0205699	3.286528	0.3105313	0.7561570	0.9919616
ASV185	22.3103569	3.5452616	3.289507	1.0777487	NA	NA
ASV186	53.4372593	-0.4119162	2.732149	-0.1507664	0.8801600	0.9919616
ASV187	18.1301853	2.8595457	3.289325	0.8693411	NA	NA
ASV188	7.9598140	-5.8981147	3.391125	-1.7392797	0.0819856	0.8011152
ASV189	53.8165630	1.1481840	2.550386	0.4502001	0.6525662	0.9919616
ASV190	70.8269217	-1.4271641	2.398874	-0.5949309	0.5518896	0.9919616
ASV191	56.6177617	-0.5870391	2.566645	-0.2287185	0.8190877	0.9919616
ASV192	46.6467565	1.0898651	2.397980	0.4544930	0.6494740	0.9919616
ASV193	48.7987065	4.9216474	3.289188	1.4963108	NA	NA
ASV194	40.7746917	1.6752866	2.573610	0.6509481	0.5150800	0.9919616
ASV195	48.1589724	1.2649562	3.055738	0.4139610	0.6789027	0.9919616
ASV196	4.3287634	6.4862999	3.325683	1.9503662	NA	NA
ASV197	54.5010051	1.2348629	3.286494	0.3757386	0.7071112	0.9919616
ASV198	50.8565292	-22.8484206	3.389657	-6.7406292	0.0000000	0.0000000
ASV199	3.8264294	6.2027020	3.326226	1.8647865	NA	NA
ASV200	47.7774942	0.3886691	3.072326	0.1265064	0.8993311	0.9919616
ASV201	44.2030046	-1.5236868	2.560120	-0.5951623	0.5517350	0.9919616
ASV202	22.6950177	0.4795746	3.288250	0.1458449	NA	NA
ASV203	54.1301506	-1.3265726	2.723668	-0.4870538	0.6262202	0.9919616
ASV204	48.8049748	1.8067306	3.286622	0.5497226	0.5825096	0.9919616
ASV205	55.2765466	-3.5158315	2.488553	-1.4128016	0.1577141	0.8239591
ASV206	26.3535301	2.2435286	3.287745	0.6823914	0.4949915	0.9919616
ASV207	38.4959885	1.0155406	3.286956	0.3089608	0.7573513	0.9919616
ASV208	14.9229529	0.4113114	3.255043	0.1263613	NA	NA
ASV209	11.5976292	2.0043164	3.290469	0.6091279	NA	NA
ASV210	18.3863950	1.2955764	3.254848	0.3980452	NA	NA
ASV211	29.9647993	0.2875438	3.025092	0.0950529	0.9242728	0.9919616
ASV212	34.0261320	-1.6277164	3.230154	-0.5039129	0.6143226	0.9919616
ASV213	46.1236192	3.1199158	3.287113	0.9491356	0.3425517	0.9919616
ASV214	15.6301085	3.3025154	3.290695	1.0035921	NA	NA
ASV215	46.9193516	-8.4568881	2.448534	-3.4538574	0.0005526	0.0112853
ASV216	34.4651476	0.3815038	3.210697	0.1188227	0.9054158	0.9919616
ASV217	53.9456636	-3.4929526	2.198408	-1.5888553	0.1120931	0.8119757
ASV218	5.5393251	-5.3755306	3.391884	-1.5848214	0.1130069	0.8119757
ASV219	35.0962495	-0.0550694	3.287560	-0.0167508	0.9866354	0.9952295
ASV220	28.2677789	0.4431015	3.050036	0.1452775	NA	NA
ASV221	52.1819645	1.1771428	3.286545	0.3581703	0.7202159	0.9919616
ASV222	22.5137246	2.2530047	3.216704	0.7004078	0.4836727	0.9919616
ASV223	48.3942311	-0.7778237	2.123283	-0.3663306	0.7141184	0.9919616
ASV224	23.3224540	-0.5419084	3.273306	-0.1655538	0.8685081	0.9919616
ASV225	17.1820961	3.2306344	3.290092	0.9819283	0.3261352	0.9919616
ASV226	13.9613476	1.4794383	3.289487	0.4497474	NA	NA
ASV227	15.4459582	0.7451223	3.253907	0.2289931	NA	NA
ASV228	18.0950799	-21.4099196	3.390150	-6.3153305	0.0000000	0.0000000
ASV229	8.9150563	1.2036123	3.291834	0.3656358	NA	NA
ASV230	18.1872632	2.2421373	3.288772	0.6817553	0.4953937	0.9919616
ASV231	36.9628833	1.2164316	3.286991	0.3700745	0.7113270	0.9919616
ASV232	8.4767189	1.2217978	3.292156	0.3711239	NA	NA
ASV233	6.3483555	-5.5716495	3.391567	-1.6427952	0.1004253	0.8011152
ASV234	14.2772380	1.1927770	3.289440	0.3626080	NA	NA
ASV235	24.2052086	0.7532608	3.013809	0.2499365	0.8026365	0.9919616
ASV236	9.2810667	24.0945809	3.323988	7.2486960	NA	NA
ASV237	36.6586500	-0.1201238	2.605021	-0.0461124	0.9632206	0.9919616
ASV238	16.3895841	2.5984305	3.289433	0.7899326	0.4295671	0.9919616
ASV239	27.5059995	1.9752325	3.207687	0.6157809	0.5380391	0.9919616
ASV240	41.7755558	-0.3194602	2.588633	-0.1234088	0.9017834	0.9919616
ASV241	37.3896336	0.5503419	2.790097	0.1972483	0.8436332	0.9919616
ASV242	9.4064437	1.0695397	3.165097	0.3379169	0.7354258	0.9919616
ASV243	21.2411784	2.5452159	3.288484	0.7739784	0.4389435	0.9919616
ASV244	49.5806522	0.0121312	3.082893	0.0039350	0.9968603	0.9968603
ASV245	10.6788339	0.6902524	3.291137	0.2097306	NA	NA
ASV246	13.0591888	0.0089209	3.290996	0.0027107	NA	NA
ASV247	9.9327054	-6.2175365	3.390779	-1.8336600	0.0667045	0.8011152
ASV248	32.5375699	3.7352951	3.288473	1.1358753	NA	NA
ASV249	35.4814086	-0.8575129	2.796340	-0.3066554	0.7591057	0.9919616
ASV250	52.1309124	0.8530368	3.286586	0.2595510	0.7952101	0.9919616
ASV251	19.5520259	2.0659194	3.288457	0.6282338	0.5298508	0.9919616
ASV252	39.2518547	-0.5221785	2.986130	-0.1748679	0.8611834	0.9919616
ASV253	43.6216423	2.6592966	3.286972	0.8090414	0.4184913	0.9919616
ASV254	25.4158867	1.2303131	3.021052	0.4072465	0.6838269	0.9919616
ASV255	12.4987439	24.4982886	3.323607	7.3709940	NA	NA
ASV256	24.2134821	-1.5551125	3.010744	-0.5165209	0.6054906	0.9919616
ASV257	39.0395996	-1.9672033	3.290153	-0.5979064	0.5499024	0.9919616
ASV258	27.0924502	0.6565059	2.833005	0.2317348	0.8167440	0.9919616
ASV259	35.6569182	-8.0610652	2.652187	-3.0394032	0.0023705	0.0459872
ASV260	16.9292489	0.0401853	3.289687	0.0122155	NA	NA
ASV261	19.3507345	1.7627943	3.288394	0.5360654	NA	NA
ASV262	34.2345214	1.4738692	3.190006	0.4620270	0.6440619	0.9919616
ASV263	29.8613336	0.6897162	3.287487	0.2098004	0.8338234	0.9919616
ASV264	32.6662047	-2.1561104	3.201494	-0.6734701	0.5006482	0.9919616
ASV265	7.6861822	23.8374323	3.324296	7.1706707	NA	NA
ASV266	14.8354575	0.7898259	3.289476	0.2401069	0.8102474	0.9919616
ASV267	31.9668768	-0.5587962	3.288257	-0.1699369	0.8650598	0.9919616
ASV268	31.9146218	-0.8801095	3.238485	-0.2717658	0.7858021	0.9919616
ASV269	11.4759410	4.5683018	3.298466	1.3849778	0.1660593	0.8239591
ASV270	20.4304992	1.6167745	3.288219	0.4916870	0.6229407	0.9919616
ASV271	4.6138978	-5.1113368	3.392386	-1.5067084	0.1318854	0.8239591
ASV272	7.5106329	1.2960818	3.292987	0.3935885	NA	NA
ASV273	20.1399236	1.9597788	3.288329	0.5959802	0.5511885	0.9919616
ASV274	25.6686758	-0.2381852	3.260728	-0.0730466	0.9417690	0.9919616
ASV275	21.9567158	1.6583223	3.288030	0.5043513	0.6140145	0.9919616
ASV276	21.5756826	-21.6536144	3.390027	-6.3874459	0.0000000	0.0000000
ASV277	34.6437893	2.3887362	2.953981	0.8086498	0.4187166	0.9919616
ASV278	45.0953610	-1.0949873	2.758066	-0.3970127	0.6913581	0.9919616
ASV279	39.5523689	-0.7212364	3.287885	-0.2193618	0.8263682	0.9919616
ASV280	31.4389091	0.7801898	3.235443	0.2411385	NA	NA
ASV281	16.6110098	1.7512755	3.288877	0.5324843	0.5943906	0.9919616
ASV282	30.2849348	0.9869484	3.287371	0.3002243	0.7640061	0.9919616
ASV283	7.4610846	1.0107439	3.293197	0.3069187	0.7589052	0.9919616
ASV284	22.4614687	0.6796489	3.256092	0.2087315	0.8346578	0.9919616
ASV285	10.2558500	-0.4690832	3.293867	-0.1424111	NA	NA
ASV286	33.6273193	-0.0593059	3.235164	-0.0183316	0.9853743	0.9952295
ASV287	10.9558645	3.4559180	3.293399	1.0493469	NA	NA
ASV288	23.1469912	-0.9516819	3.252442	-0.2926054	0.7698238	0.9919616
ASV289	33.7411213	0.5850212	3.287291	0.1779645	0.8587509	0.9919616
ASV290	5.0439992	-5.2372063	3.392136	-1.5439260	0.1226063	0.8239591
ASV291	17.7192392	0.2909657	3.289209	0.0884607	0.9295105	0.9919616
ASV292	16.3203729	1.6685669	3.288923	0.5073293	0.6119238	0.9919616
ASV293	21.7574716	-0.5882902	3.263430	-0.1802674	0.8569426	0.9919616
ASV294	16.4221681	2.7611878	3.289603	0.8393681	0.4012628	0.9919616
ASV295	6.0180935	-5.4934464	3.391688	-1.6196791	0.1053012	0.8011152
ASV296	7.8064369	3.9845315	3.299557	1.2075958	NA	NA
ASV297	9.9510813	24.1891958	3.323889	7.2773777	NA	NA
ASV298	1.4781120	4.9375253	3.331789	1.4819440	NA	NA
ASV299	14.9884137	-21.1466952	3.369857	-6.2752507	0.0000000	0.0000000
ASV300	31.0897454	1.4297142	3.287267	0.4349249	0.6636170	0.9919616
ASV301	7.9523599	23.8837630	3.324237	7.1847365	NA	NA
ASV302	18.2433368	-0.8310908	3.290990	-0.2525352	0.8006274	0.9919616
ASV303	2.0995546	-3.9849366	3.395931	-1.1734444	0.2406176	0.8845429
ASV304	23.9470950	1.2435789	3.287824	0.3782377	NA	NA
ASV305	19.4090722	2.0256975	3.288461	0.6160018	0.5378934	0.9919616
ASV306	6.8574877	1.2261718	3.293736	0.3722739	NA	NA
ASV307	25.5499679	1.5054092	3.025120	0.4976362	0.6187405	0.9919616
ASV308	15.3383733	0.1406839	2.996615	0.0469476	0.9625550	0.9919616
ASV309	20.5325467	-4.2107017	2.789221	-1.5096336	0.1311369	0.8239591
ASV310	12.6632098	-0.2869019	2.970557	-0.0965818	NA	NA
ASV311	6.2437627	-5.5469686	3.391604	-1.6354999	0.1019443	0.8011152
ASV312	10.0826985	1.8543329	3.291134	0.5634328	NA	NA
ASV313	22.3408421	-0.0531274	3.288762	-0.0161542	NA	NA
ASV314	7.8861252	2.6442201	3.293810	0.8027847	NA	NA
ASV315	15.1683373	1.9318736	3.289260	0.5873277	NA	NA
ASV316	13.3758267	3.6666508	3.292579	1.1136104	0.2654464	0.9179515
ASV317	8.4905960	1.2789185	3.292129	0.3884777	NA	NA
ASV318	30.7943441	0.6231298	3.239203	0.1923714	0.8474513	0.9919616
ASV319	28.2112952	-0.2694610	2.848642	-0.0945928	NA	NA
ASV320	18.4037310	2.0780666	3.288649	0.6318907	NA	NA
ASV321	22.5055690	2.6995797	3.288432	0.8209323	0.4116848	0.9919616
ASV322	14.8047208	-2.2036178	3.244258	-0.6792363	0.4969881	0.9919616
ASV323	9.4589579	0.8903929	3.291664	0.2704993	NA	NA
ASV324	21.0513109	1.6556005	3.229021	0.5127252	0.6081435	0.9919616
ASV325	6.4577070	1.6873393	3.294222	0.5122118	NA	NA
ASV326	22.5854420	0.0478903	3.252817	0.0147227	0.9882534	0.9952295
ASV327	19.9823089	-2.8710827	2.984753	-0.9619164	0.3360916	0.9919616
ASV328	19.2104796	-1.3709899	2.624961	-0.5222896	0.6014687	0.9919616
ASV329	5.0530925	6.7096768	3.325228	2.0178098	0.0436111	0.7356998
ASV330	25.7234206	-21.8979274	3.388500	-6.4624244	0.0000000	0.0000000
ASV331	7.2275408	4.6623864	3.307107	1.4098081	0.1585964	0.8239591
ASV332	15.4107524	1.7090247	3.289137	0.5195966	0.6033448	0.9919616
ASV333	9.8369510	2.2011020	3.291550	0.6687129	NA	NA
ASV334	31.8099597	-22.1947664	3.351546	-6.6222471	0.0000000	0.0000000
ASV335	27.1051847	2.0826697	3.229479	0.6448934	0.5189963	0.9919616
ASV336	6.5372984	-5.6122495	3.391506	-1.6547956	0.0979660	0.8011152
ASV337	2.5298948	-4.2431824	3.394860	-1.2498844	0.2113418	0.8814503
ASV338	9.6617263	3.5721724	3.294911	1.0841482	0.2782991	0.9229064
ASV339	25.5256802	0.0969723	3.288212	0.0294909	0.9764731	0.9919616
ASV340	3.2319751	6.0645677	3.326761	1.8229647	0.0683087	0.8011152
ASV341	22.3250017	1.1194853	3.256602	0.3437587	0.7310277	0.9919616
ASV342	6.8223288	0.7725510	3.294216	0.2345174	NA	NA
ASV343	2.2755830	5.5601056	3.328540	1.6704339	NA	NA
ASV344	11.7327809	0.5186811	3.290811	0.1576150	0.8747602	0.9919616
ASV345	10.3763412	24.2051116	3.323832	7.2822900	NA	NA
ASV346	13.5946642	0.5370045	3.290062	0.1632202	0.8703451	0.9919616
ASV347	3.5719895	22.7924678	3.326353	6.8520890	NA	NA
ASV348	10.9576087	1.1924889	3.290646	0.3623875	0.7170625	0.9919616
ASV349	12.6577918	2.3155956	3.290286	0.7037673	0.4815777	0.9919616
ASV350	25.2809479	0.5876692	3.261714	0.1801719	0.8570176	0.9919616
ASV351	1.0759840	-3.0155157	3.402181	-0.8863478	0.3754301	0.9919616
ASV352	2.5994486	-4.2911586	3.394681	-1.2640830	0.2062002	0.8791834
ASV353	8.4659340	-0.2966824	3.294978	-0.0900408	NA	NA
ASV354	6.6175349	-5.6308267	3.391479	-1.6602865	0.0968568	0.8011152
ASV355	9.6891171	0.5318415	3.291923	0.1615595	NA	NA
ASV356	16.1466007	2.4615437	3.289362	0.7483348	NA	NA
ASV357	1.1744862	1.4962560	3.333397	0.4488682	0.6535267	0.9919616
ASV358	4.5719080	-5.0969026	3.392416	-1.5024402	0.1329835	0.8239591
ASV359	4.1099175	-4.9429479	3.392757	-1.4569117	0.1451408	0.8239591
ASV360	9.1161164	0.2616687	3.292802	0.0794669	NA	NA
ASV361	2.9150425	5.9160224	3.327222	1.7780667	NA	NA
ASV362	10.5798567	2.2664968	3.291188	0.6886562	0.4910397	0.9919616
ASV363	8.1041759	3.5559606	2.750697	1.2927492	0.1960978	0.8646131
ASV364	4.3866723	-5.0390297	3.392540	-1.4853264	0.1374574	0.8239591
ASV365	17.5116212	2.8908571	3.289501	0.8788132	NA	NA
ASV366	11.1577084	-0.1605420	3.292323	-0.0487625	NA	NA
ASV367	15.7930769	0.1635221	2.781788	0.0587831	0.9531249	0.9919616
ASV368	22.6282192	2.7141937	3.288428	0.8253773	0.4091574	0.9919616
ASV369	5.7247150	1.3677967	3.295281	0.4150774	NA	NA
ASV370	1.1381014	-3.0806698	3.401618	-0.9056484	0.3651220	0.9919616
ASV371	5.1695403	1.9520955	3.296663	0.5921428	NA	NA
ASV372	4.4594402	23.0867014	3.325592	6.9421330	NA	NA
ASV373	2.1918493	-4.0429070	3.395674	-1.1906053	0.2338086	0.8845429
ASV374	15.9080229	1.6743958	3.289014	0.5090875	0.6106909	0.9919616
ASV375	2.3550960	-4.1492944	3.395227	-1.2220962	0.2216713	0.8845429
ASV376	6.9135138	-5.6946465	3.391388	-1.6791490	0.0931230	0.8011152
ASV377	8.1874175	1.4892460	3.292338	0.4523369	NA	NA
ASV378	2.4720810	5.6770332	3.328070	1.7058032	NA	NA
ASV379	9.7603254	3.7740186	3.295689	1.1451381	NA	NA
ASV380	3.7848130	6.2923456	3.326140	1.8917862	NA	NA
ASV381	8.7680420	2.1490889	3.292217	0.6527786	NA	NA
ASV382	9.7940068	3.4540112	3.294331	1.0484712	NA	NA
ASV383	23.3290658	-7.4488281	3.143582	-2.3695349	0.0178105	0.3290697
ASV384	8.8949448	1.6494805	3.291812	0.5010859	NA	NA
ASV385	3.1147029	6.0115655	3.326920	1.8069463	NA	NA
ASV386	6.6812132	-5.6445462	3.391459	-1.6643415	0.0960442	0.8011152
ASV387	6.4190207	0.1818046	3.296135	0.0551569	NA	NA
ASV388	6.4079418	-5.5823730	3.391551	-1.6459649	0.0997710	0.8011152
ASV389	11.3484217	2.7554524	3.291442	0.8371566	0.4025046	0.9919616
ASV390	7.0834173	1.6255554	3.293411	0.4935780	0.6216042	0.9919616
ASV391	6.8358996	-5.6774673	3.391412	-1.6740716	0.0941165	0.8011152
ASV392	3.1347386	-4.5518680	3.393806	-1.3412279	0.1798465	0.8239591
ASV393	10.2962006	-6.2675505	3.390732	-1.8484359	0.0645393	0.8011152
ASV394	3.4740917	6.1690542	3.326464	1.8545382	NA	NA
ASV395	10.2387843	0.5422407	3.291565	0.1647364	NA	NA
ASV396	23.2932545	1.2940604	3.287885	0.3935845	0.6938879	0.9919616
ASV397	4.3745643	2.4271239	3.299739	0.7355502	NA	NA
ASV398	4.7484277	5.4350098	3.324544	1.6348136	NA	NA
ASV399	3.5911397	-4.7489913	3.393242	-1.3995438	0.1616500	0.8239591
ASV400	13.8403938	1.3775190	3.289532	0.4187583	0.6753928	0.9919616
ASV401	7.1472598	1.7943223	3.293413	0.5448216	NA	NA
ASV402	6.0998737	-5.5128067	3.391657	-1.6254020	0.1040769	0.8011152
ASV403	7.1915050	-1.5644982	3.167674	-0.4938950	NA	NA
ASV404	14.0543729	0.9313149	3.289614	0.2831077	0.7770943	0.9919616
ASV405	13.8193594	-6.6931086	3.358081	-1.9931350	0.0462467	0.7476544
ASV406	3.3270091	-1.4441478	3.326308	-0.4341594	0.6641727	0.9919616
ASV407	9.7024443	1.6426026	3.291280	0.4990771	0.6177251	0.9919616
ASV408	11.7020946	-0.6230667	3.293328	-0.1891906	0.8499434	0.9919616
ASV409	8.9737252	1.9195733	3.291875	0.5831247	0.5598094	0.9919616
ASV410	2.7476842	5.8320359	3.327504	1.7526758	NA	NA
ASV411	9.7297698	0.1064829	3.292665	0.0323394	0.9742014	0.9919616
ASV412	6.9994500	-5.7105018	3.391367	-1.6838351	0.0922135	0.8011152
ASV413	0.0000000	0.0000000	0.000000	0.0000000	1.0000000	NA
ASV414	4.4228487	1.0500917	3.298443	0.3183598	NA	NA
ASV415	2.5690254	5.7325532	3.327861	1.7225941	NA	NA
ASV416	12.3453881	-0.6414907	3.292975	-0.1948058	0.8455449	0.9919616
ASV417	6.6150192	-5.6323120	3.391477	-1.6607255	0.0967686	0.8011152
ASV418	3.8036444	-4.8298986	3.393032	-1.4234757	0.1545983	0.8239591
ASV419	16.0023439	2.6099305	3.289538	0.7934033	0.4275429	0.9919616
ASV420	2.4239814	5.6495918	3.328177	1.6975033	NA	NA
ASV421	2.3881343	-4.1661621	3.395160	-1.2270887	0.2197892	0.8845429
ASV422	6.7740307	-5.6656290	3.391429	-1.6705727	0.0948061	0.8011152
ASV423	7.6803861	-5.8441031	3.391191	-1.7233186	0.0848309	0.8011152
ASV424	4.3557170	0.9889644	3.298775	0.2997975	0.7643316	0.9919616
ASV425	5.7186088	6.8864402	3.324914	2.0711635	0.0383435	0.6762403
ASV426	8.5202847	3.6565795	3.296575	1.1092055	0.2673415	0.9179515
ASV427	2.5642318	0.3011807	3.311220	0.0909576	NA	NA
ASV428	8.2709097	-5.9519622	3.391061	-1.7551916	0.0792265	0.8011152
ASV429	5.3842392	2.7507337	3.298058	0.8340464	NA	NA
ASV430	3.3595531	6.1204137	3.326599	1.8398409	NA	NA
ASV431	13.5419296	2.0505148	3.289780	0.6232985	0.5330884	0.9919616
ASV432	2.1745387	-4.0352213	3.395707	-1.1883301	0.2347034	0.8845429
ASV433	3.4265072	-4.6755238	3.393444	-1.3778110	0.1682617	0.8239591
ASV434	2.1011621	-3.9762197	3.395971	-1.1708639	0.2416535	0.8845429
ASV435	1.9410121	-3.8718021	3.396464	-1.1399508	0.2543068	0.9136207
ASV436	6.1651223	4.6268684	3.310374	1.3976874	NA	NA
ASV437	7.9489448	0.7476275	3.293017	0.2270342	0.8203971	0.9919616
ASV438	3.9901374	2.2086291	3.300447	0.6691909	NA	NA
ASV439	1.4523470	-3.4419637	3.398915	-1.0126654	0.3112200	0.9759373
ASV440	1.1979627	4.6338383	3.333960	1.3898904	NA	NA
ASV441	3.0751319	-4.5297205	3.393875	-1.3346752	0.1819827	0.8239591
ASV442	7.7704493	-0.5019681	3.296708	-0.1522634	0.8789792	0.9919616
ASV443	3.1325607	6.0167761	3.326904	1.8085211	NA	NA
ASV444	7.8490268	1.7870771	3.292708	0.5427377	NA	NA
ASV445	0.4196783	-1.7327632	3.419279	-0.5067628	0.6123213	0.9919616
ASV446	11.1654858	-6.3848205	3.390627	-1.8830798	0.0596896	0.8011152
ASV447	6.7424884	4.1324456	3.303034	1.2511060	NA	NA
ASV448	3.3273285	-4.6401090	3.393544	-1.3673341	0.1715206	0.8239591
ASV449	3.7199159	6.2651016	3.326209	1.8835562	NA	NA
ASV450	4.5030560	6.5411053	3.325565	1.9669157	NA	NA
ASV451	2.2228125	-4.0566065	3.395614	-1.1946605	0.2322197	0.8845429
ASV452	7.1174911	1.0877348	3.293514	0.3302657	NA	NA
ASV453	8.1716461	2.9002726	3.294150	0.8804314	0.3786257	0.9919616
ASV454	10.3049731	0.8196009	3.291217	0.2490267	NA	NA
ASV455	7.3939638	3.2358032	3.296264	0.9816578	NA	NA
ASV456	0.7229418	3.9046857	3.341463	1.1685556	NA	NA
ASV457	5.3749674	2.4530743	3.297154	0.7439975	NA	NA
ASV458	3.9974429	22.8884702	3.325947	6.8817911	NA	NA
ASV459	2.6255649	5.7665333	3.327736	1.7328698	NA	NA
ASV460	3.8859625	-4.8643999	3.392946	-1.4336803	0.1516635	0.8239591
ASV461	4.3263698	-5.0178655	3.392587	-1.4790676	0.1391222	0.8239591
ASV462	2.2113514	5.5171646	3.328722	1.6574424	NA	NA
ASV463	5.6843428	-5.4116591	3.391822	-1.5955019	0.1106000	0.8119757
ASV464	15.0021654	1.3422671	3.289220	0.4080806	0.6832145	0.9919616
ASV465	7.9773806	1.2336232	3.292576	0.3746681	NA	NA
ASV466	1.8739961	-3.8099752	3.396774	-1.1216453	0.2620133	0.9179515
ASV467	11.8041325	2.4503960	3.290782	0.7446244	0.4564988	0.9919616
ASV468	4.0685959	-4.9270168	3.392795	-1.4522000	0.1464460	0.8239591
ASV469	8.8361550	-6.0470282	3.390954	-1.7832821	0.0745404	0.8011152
ASV470	3.2745011	6.0834242	3.326706	1.8286632	NA	NA
ASV471	12.8718050	1.8700006	3.289916	0.5684038	0.5697608	0.9919616
ASV472	5.1334321	-5.2625545	3.392088	-1.5514205	0.1208009	0.8239591
ASV473	4.7372088	3.9529295	3.308363	1.1948294	NA	NA
ASV474	1.3573802	1.4857576	3.326681	0.4466186	NA	NA
ASV475	0.2235279	-0.9057717	3.424764	-0.2644771	0.7914123	0.9919616
ASV476	3.0699998	-4.5230462	3.393896	-1.3327005	0.1826301	0.8239591
ASV477	2.4847749	-4.2228624	3.394937	-1.2438705	0.2135472	0.8814503
ASV478	10.2826705	1.3143547	3.290955	0.3993840	NA	NA
ASV479	6.6547053	1.9096359	3.294104	0.5797134	NA	NA
ASV480	10.4738435	1.7826559	3.290896	0.5416932	NA	NA
ASV481	0.2249523	-1.1882658	3.423912	-0.3470492	0.7285544	0.9919616
ASV482	2.5087220	-4.2308116	3.394907	-1.2462232	0.2126825	0.8814503
ASV483	3.6658572	-4.7790901	3.393163	-1.4084471	0.1589987	0.8239591
ASV484	11.2542792	2.5356928	3.291159	0.7704558	0.4410296	0.9919616
ASV485	6.1345981	23.5259778	3.324751	7.0760117	NA	NA
ASV486	9.2624884	-6.1148235	3.390882	-1.8033134	0.0713390	0.8011152
ASV487	0.8368765	-2.6543018	3.405805	-0.7793463	0.4357757	0.9919616
ASV488	4.9177500	-5.1987673	3.392210	-1.5325607	0.1253841	0.8239591
ASV489	8.1782436	0.2400018	3.293691	0.0728671	0.9419118	0.9919616
ASV490	1.7402000	5.1741333	3.330387	1.5536130	NA	NA
ASV491	5.8695864	1.1652273	3.295163	0.3536175	NA	NA
ASV492	4.7635678	4.0566959	3.309351	1.2258284	NA	NA
ASV493	12.1468607	0.9961102	3.290228	0.3027481	0.7620818	0.9919616
ASV494	3.3296712	6.1075470	3.326636	1.8359529	0.0663646	0.8011152
ASV495	5.5991049	3.3507677	3.300292	1.0152943	NA	NA
ASV496	8.3510393	2.2072796	3.292637	0.6703683	0.5026230	0.9919616
ASV497	11.3134399	5.4331437	3.307194	1.6428258	0.1004190	0.8011152
ASV498	3.3568899	-4.6533701	3.393506	-1.3712572	0.1702948	0.8239591
ASV499	7.8805435	2.4221323	3.293401	0.7354501	NA	NA
ASV500	2.0580684	3.6708641	3.327078	1.1033296	NA	NA
ASV501	0.1108584	1.3643242	3.418344	0.3991184	0.6898059	0.9919616
ASV502	3.3367237	2.3049513	3.303658	0.6976967	NA	NA
ASV503	2.3020737	-4.1059955	3.395405	-1.2092800	0.2265553	0.8845429
ASV504	2.1998060	2.3223581	3.313194	0.7009424	NA	NA
ASV505	1.5662804	5.0142549	3.331309	1.5051904	NA	NA
ASV506	4.6581141	0.8176914	3.298194	0.2479210	NA	NA
ASV507	2.8067152	5.8610626	3.327405	1.7614517	NA	NA
ASV508	1.8814542	5.2848882	3.329806	1.5871460	NA	NA
ASV509	4.7639619	-5.1581567	3.392290	-1.5205528	0.1283721	0.8239591
ASV510	0.3779434	-1.5609347	3.421622	-0.4561973	0.6482481	0.9919616
ASV511	0.6786105	3.8101904	3.342741	1.1398400	NA	NA
ASV512	3.4953751	-4.7085942	3.393352	-1.3875939	0.1652607	0.8239591
ASV513	2.1844963	1.0296233	3.311992	0.3108773	NA	NA
ASV514	2.2897368	5.5693694	3.328501	1.6732365	NA	NA
ASV515	3.2428620	-4.6014779	3.393657	-1.3559053	0.1751293	0.8239591
ASV516	0.7765836	0.4665551	3.365866	0.1386137	0.8897554	0.9919616
ASV517	3.1752097	0.0522770	3.307868	0.0158038	NA	NA
ASV518	9.8187479	3.9116049	3.296317	1.1866591	0.2353621	0.8845429
ASV519	1.3976231	4.8560517	3.332328	1.4572552	NA	NA
ASV520	2.4468189	1.8987541	3.308754	0.5738578	NA	NA
ASV521	1.4049190	1.5595255	3.324962	0.4690357	NA	NA
ASV522	6.9852153	3.3538173	3.217753	1.0422854	0.2972794	0.9568046
ASV523	0.8317690	-2.6343256	3.406033	-0.7734292	0.4392684	0.9919616
ASV524	2.4756094	1.0398993	3.308792	0.3142837	NA	NA
ASV525	2.3576093	-4.1442696	3.395248	-1.2206089	0.2222341	0.8845429
ASV526	2.6862027	-4.3273613	3.394550	-1.2747968	0.2023812	0.8791834
ASV527	5.4731279	6.4792147	3.325094	1.9485809	NA	NA
ASV528	1.9096601	-3.8351181	3.396646	-1.1290897	0.2588600	0.9179515
ASV529	4.7518282	1.8717140	3.297526	0.5676116	NA	NA
ASV530	2.6115207	2.4179023	3.309369	0.7306234	NA	NA
ASV531	0.7152892	-2.4219293	3.408660	-0.7105224	0.4773802	0.9919616
ASV532	3.8696715	3.2796212	3.306362	0.9919124	NA	NA
ASV533	0.9355717	4.2764928	3.337169	1.2814732	NA	NA
ASV534	3.6924976	4.2743293	3.319008	1.2878335	NA	NA
ASV535	3.1687303	1.3049795	3.303241	0.3950603	NA	NA
ASV536	0.6086405	-2.2016835	3.411820	-0.6453106	0.5187259	0.9919616
ASV537	5.8014094	1.2609867	3.295224	0.3826710	0.7019637	0.9919616
ASV538	1.9427797	-0.4474780	3.329109	-0.1344137	0.8930754	0.9919616
ASV539	3.9430191	-4.8848708	3.392896	-1.4397350	0.1499424	0.8239591
ASV540	1.2851406	1.8692532	3.329799	0.5613712	NA	NA
ASV541	4.4770362	5.5667952	3.325004	1.6742220	0.0940870	0.8011152
ASV542	0.6891409	-2.3415078	3.409759	-0.6867077	0.4922669	0.9919616
ASV543	0.8174920	-2.6208135	3.406189	-0.7694269	0.4416399	0.9919616
ASV544	3.2208127	-4.5886563	3.393695	-1.3521121	0.1763395	0.8239591
ASV545	1.1798472	-3.1420619	3.401110	-0.9238343	0.3555726	0.9919616
ASV546	1.2847054	1.3156473	3.329223	0.3951815	NA	NA
ASV547	2.2870624	3.1179216	3.318576	0.9395361	0.3474556	0.9919616
ASV548	1.8426105	-3.7853417	3.396901	-1.1143517	0.2651284	0.9179515
ASV549	4.1550795	2.6360726	3.301294	0.7984969	0.4245822	0.9919616
ASV550	1.1482260	2.2035884	3.335215	0.6607036	NA	NA
ASV551	2.6363579	-4.3017006	3.394642	-1.2672028	0.2050828	0.8791834
ASV552	0.9606208	0.9642660	3.346299	0.2881590	0.7732250	0.9919616
ASV553	0.4904125	-1.8835702	3.416901	-0.5512510	0.5814617	0.9919616
ASV554	1.6106202	0.5682425	3.323866	0.1709583	NA	NA
ASV555	0.4349425	1.8073549	3.379080	0.5348661	0.5927425	0.9919616
ASV556	0.6202405	-0.1226365	3.397837	-0.0360925	0.9712086	0.9919616
ASV557	0.8806597	2.1304638	3.343906	0.6371183	0.5240478	0.9919616

#umbral
alpha = 0.05
# Ordenamos la tabla de resultados
res.PR.sCR  = res.PR.sCR [order(res.PR.sCR $padj, na.last=NA), ]
# Filtramos según nuestro umbral alpha
sigtab = res.PR.sCR [(res.PR.sCR $padj < alpha), ]
# Le agregamos la taxonomía a la tabla
sigtab = cbind(as(sigtab, "data.frame"), as(tax_table(psd5)[rownames(sigtab), ], "matrix"))

# Manipulaciones varias para finalmente graficar los resultados
sigtabgen = subset(sigtab, !is.na(Genus))
# Phylum order
x = tapply(sigtabgen$log2FoldChange, sigtabgen$Phylum, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Phylum = factor(as.character(sigtabgen$Phylum), levels=names(x))
# Genus order
x = tapply(sigtabgen$log2FoldChange, sigtabgen$Genus, function(x) max(x))
x = sort(x, TRUE)
sigtabgen$Genus = factor(as.character(sigtabgen$Genus), levels=names(x))
ggplot(sigtabgen, aes(y=Genus, x=log2FoldChange, color=Phylum)) + 
    geom_vline(xintercept = 0.0, color = "gray", size = 0.5) +
    geom_point(size=4) + 
    theme(axis.text.x = element_text(angle = -90, hjust = 0, vjust=0.5, size = 10), axis.text.y = element_text(size = 13), legend.text = element_text(size = 13) )

df_res <- cbind(as(res.PR.sCR, "data.frame"), as(tax_table(psd5)[rownames(res.PR.sCR), ], "matrix"))
df_res <- df_res %>%
  rownames_to_column(var = "OTU") %>%
  arrange(padj)

(fdr_otu <- df_res %>%
  dplyr::filter(padj < 0.05))

##       OTU  baseMean log2FoldChange    lfcSE      stat       pvalue         padj
## 1  ASV143  90.57919     -23.647542 2.798983 -8.448620 2.947668e-17 1.143695e-14
## 2   ASV78 138.27535     -24.242780 3.338224 -7.262179 3.809041e-13 3.986267e-11
## 3  ASV105 112.93822     -23.953855 3.303113 -7.251903 4.109554e-13 3.986267e-11
## 4  ASV110 108.77143     -23.903707 3.290335 -7.264825 3.735217e-13 3.986267e-11
## 5  ASV198  50.85653     -22.848421 3.389657 -6.740629 1.577022e-11 1.223769e-09
## 6  ASV334  31.80996     -22.194766 3.351546 -6.622247 3.537789e-11 2.287770e-09
## 7  ASV330  25.72342     -21.897927 3.388500 -6.462424 1.030386e-10 5.711284e-09
## 8   ASV48  21.06044     -21.618912 3.390043 -6.377180 1.803782e-10 7.063576e-09
## 9   ASV61  20.96650     -21.614137 3.390046 -6.375766 1.820509e-10 7.063576e-09
## 10 ASV276  21.57568     -21.653614 3.390027 -6.387446 1.686793e-10 7.063576e-09
## 11  ASV95  18.96324     -21.466058 3.390115 -6.331955 2.420742e-10 8.046069e-09
## 12 ASV173  27.80827     -21.431890 3.389883 -6.322310 2.576820e-10 8.046069e-09
## 13 ASV228  18.09508     -21.409920 3.390150 -6.315330 2.695848e-10 8.046069e-09
## 14 ASV299  14.98841     -21.146695 3.369857 -6.275251 3.490714e-10 9.674264e-09
## 15 ASV128  14.00293     -21.052789 3.390374 -6.209577 5.312756e-10 1.374233e-08
## 16 ASV165  26.28655     -20.722081 3.389912 -6.112868 9.785659e-10 2.373022e-08
## 17 ASV147  76.99358      -9.171426 2.113207 -4.340051 1.424494e-05 3.251197e-04
## 18 ASV171  65.75344      -8.943829 2.569298 -3.481040 4.994711e-04 1.076638e-02
## 19 ASV215  46.91935      -8.456888 2.448534 -3.453857 5.526295e-04 1.128528e-02
## 20 ASV259  35.65692      -8.061065 2.652187 -3.039403 2.370474e-03 4.598719e-02
##     Kingdom           Phylum          Class             Order
## 1  Bacteria Actinobacteriota Actinobacteria Bifidobacteriales
## 2  Bacteria Actinobacteriota Actinobacteria Bifidobacteriales
## 3  Bacteria Actinobacteriota Actinobacteria Bifidobacteriales
## 4  Bacteria Actinobacteriota Actinobacteria Bifidobacteriales
## 5  Bacteria     Bacteroidota    Bacteroidia     Bacteroidales
## 6  Bacteria       Firmicutes     Clostridia    Lachnospirales
## 7  Bacteria Actinobacteriota Coriobacteriia  Coriobacteriales
## 8  Bacteria       Firmicutes        Bacilli   Lactobacillales
## 9  Bacteria       Firmicutes        Bacilli   Lactobacillales
## 10 Bacteria       Firmicutes     Clostridia   Oscillospirales
## 11 Bacteria       Firmicutes        Bacilli   Lactobacillales
## 12 Bacteria     Bacteroidota    Bacteroidia     Bacteroidales
## 13 Bacteria     Bacteroidota    Bacteroidia     Bacteroidales
## 14 Bacteria       Firmicutes     Clostridia   Oscillospirales
## 15 Bacteria       Firmicutes        Bacilli   Lactobacillales
## 16 Bacteria     Bacteroidota    Bacteroidia     Bacteroidales
## 17 Bacteria       Firmicutes        Bacilli   Lactobacillales
## 18 Bacteria       Firmicutes        Bacilli   Lactobacillales
## 19 Bacteria       Firmicutes        Bacilli   Lactobacillales
## 20 Bacteria       Firmicutes     Clostridia    Lachnospirales
##                Family                        Genus
## 1  Bifidobacteriaceae              Bifidobacterium
## 2  Bifidobacteriaceae              Bifidobacterium
## 3  Bifidobacteriaceae              Bifidobacterium
## 4  Bifidobacteriaceae              Bifidobacterium
## 5      Bacteroidaceae                  Bacteroides
## 6     Lachnospiraceae                  Coprococcus
## 7     Eggerthellaceae             Senegalimassilia
## 8    Lactobacillaceae            Ligilactobacillus
## 9    Lactobacillaceae            Ligilactobacillus
## 10    Ruminococcaceae             Faecalibacterium
## 11   Lactobacillaceae            Ligilactobacillus
## 12     Bacteroidaceae                  Bacteroides
## 13     Bacteroidaceae                  Bacteroides
## 14    Ruminococcaceae             Faecalibacterium
## 15   Lactobacillaceae            Ligilactobacillus
## 16     Bacteroidaceae                  Bacteroides
## 17   Streptococcaceae                Streptococcus
## 18   Streptococcaceae                Streptococcus
## 19   Streptococcaceae                Streptococcus
## 20    Lachnospiraceae [Ruminococcus] torques group

#Metadata counts
as(sample_data(psd5), "data.frame") -> aaa
library(plyr)
counts <- ddply(aaa, .(aaa$Respuesta_100, aaa$Sexo), nrow)
names(counts) <- c("Respuesta_100", "Sexo", "Freq")
counts

##   Respuesta_100   Sexo Freq
## 1            PR HOMBRE    4
## 2            PR  MUJER    2
## 3           sCR HOMBRE   12
## 4           sCR  MUJER    5

list(fdr_otu$OTU)

## [[1]]
##  [1] "ASV143" "ASV78"  "ASV105" "ASV110" "ASV198" "ASV334" "ASV330" "ASV48" 
##  [9] "ASV61"  "ASV276" "ASV95"  "ASV173" "ASV228" "ASV299" "ASV128" "ASV165"
## [17] "ASV147" "ASV171" "ASV215" "ASV259"

#Create a phyloseq object OTUS significatives



GP.chl = prune_taxa(taxa_names(psd5) %in% c("ASV143", "ASV78",  "ASV105", "ASV110", "ASV198", "ASV334", "ASV330", "ASV48",  "ASV61","ASV276", "ASV95",  "ASV173", "ASV228", "ASV299", "ASV128", "ASV165", "ASV147", "ASV171", "ASV215", "ASV259"), psd5)

GP.chl

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 20 taxa and 23 samples ]
## sample_data() Sample Data:       [ 23 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 20 taxa by 6 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 20 tips and 19 internal nodes ]

taxa_sums(GP.chl)

##  ASV48  ASV61  ASV78  ASV95 ASV105 ASV110 ASV128 ASV143 ASV147 ASV165 ASV171 
##   3673   3338   2556   2609   2089   2027   1890   1698   1729   1885   1575 
## ASV173 ASV198 ASV215 ASV228 ASV259 ASV276 ASV299 ASV330 ASV334 
##   1849   1483   1072   1176    861    436    355    686    647

#Comparation of abundance
taxa_names(GP.chl)

##  [1] "ASV48"  "ASV61"  "ASV78"  "ASV95"  "ASV105" "ASV110" "ASV128" "ASV143"
##  [9] "ASV147" "ASV165" "ASV171" "ASV173" "ASV198" "ASV215" "ASV228" "ASV259"
## [17] "ASV276" "ASV299" "ASV330" "ASV334"

#metadata, taxonomíaa y abundancia del taxon más abundante de cada muestra

df2 <- psmelt(GP.chl)

head(df2) %>%
  kable(format = "html", col.names = colnames(df2)) %>%
  kable_styling() %>%
  kableExtra::scroll_box(width = "100%", height = "400px")

	OTU	Sample	Abundance	SampleID	Edad_a_la_inclusion	Sexo	Respuesta_100	aloTPH_previo_si_no	Citogenética_alto_riesgo	Afect_Extramedular_PET_screening_si_no	Dias_producción	CRS_si_no	Grado_máx_CRS	HOSPITAL	Visit	Kingdom	Phylum	Class	Order	Family	Genus
379	ASV48	HCB12-V2	1843	HCB12-V2	64	HOMBRE	sCR	NO	SI	NO	110	NO	NA	HCB	V2	Bacteria	Firmicutes	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus
401	ASV61	HCB12-V2	1778	HCB12-V2	64	HOMBRE	sCR	NO	SI	NO	110	NO	NA	HCB	V2	Bacteria	Firmicutes	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus
381	ASV48	HVA04-V2	1272	HVA04-V2	74	MUJER	sCR	NO	NO	NO	100	SI	10	HVA	V2	Bacteria	Firmicutes	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus
442	ASV95	HCB12-V2	1115	HCB12-V2	64	HOMBRE	sCR	NO	SI	NO	110	NO	NA	HCB	V2	Bacteria	Firmicutes	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus
134	ASV165	ISBS03-V2	1061	ISBS03-V2	70	MUJER	sCR	NO		SI	100	SI	10	ISBS	V2	Bacteria	Bacteroidota	Bacteroidia	Bacteroidales	Bacteroidaceae	Bacteroides
399	ASV61	HVA04-V2	1029	HVA04-V2	74	MUJER	sCR	NO	NO	NO	100	SI	10	HVA	V2	Bacteria	Firmicutes	Bacilli	Lactobacillales	Lactobacillaceae	Ligilactobacillus

#Filtrar OTUS sifnificativas
xx <- df %>%
  group_by(OTU)%>%
  filter(OTU %in% c("ASV143", "ASV78",  "ASV105", "ASV110", "ASV198", "ASV334", "ASV330", "ASV48",  "ASV61","ASV276", "ASV95",  "ASV173", "ASV228", "ASV299", "ASV128", "ASV165", "ASV147", "ASV171", "ASV215", "ASV259"))

xx

## # A tibble: 460 × 21
## # Groups:   OTU [20]
##    OTU    Sample    Abundance SampleID  Edad_a_la_inclusion Sexo   Respuesta_100
##    <chr>  <chr>         <int> <chr>                   <int> <chr>  <fct>        
##  1 ASV48  HCB12-V2       1843 HCB12-V2                   64 HOMBRE sCR          
##  2 ASV61  HCB12-V2       1778 HCB12-V2                   64 HOMBRE sCR          
##  3 ASV48  HVA04-V2       1272 HVA04-V2                   74 MUJER  sCR          
##  4 ASV95  HCB12-V2       1115 HCB12-V2                   64 HOMBRE sCR          
##  5 ASV165 ISBS03-V2      1061 ISBS03-V2                  70 MUJER  sCR          
##  6 ASV61  HVA04-V2       1029 HVA04-V2                   74 MUJER  sCR          
##  7 ASV95  HVA04-V2       1020 HVA04-V2                   74 MUJER  sCR          
##  8 ASV173 ISBS03-V2      1003 ISBS03-V2                  70 MUJER  sCR          
##  9 ASV198 ISBS03-V2       819 ISBS03-V2                  70 MUJER  sCR          
## 10 ASV128 HCB12-V2        780 HCB12-V2                   64 HOMBRE sCR          
## # … with 450 more rows, and 14 more variables: aloTPH_previo_si_no <chr>,
## #   Citogenética_alto_riesgo <chr>,
## #   Afect_Extramedular_PET_screening_si_no <chr>, Dias_producción <int>,
## #   CRS_si_no <chr>, Grado_máx_CRS <int>, HOSPITAL <chr>, Visit <chr>,
## #   Kingdom <chr>, Phylum <chr>, Class <chr>, Order <chr>, Family <chr>,
## #   Genus <chr>

#Save table in .csv file
write_csv(xx, file="xx.csv")

1_Phyloseq

Eduardo

15/1/2022

Table of read counts